# Operator Spreading and the Absence of Many-Body Localization

A. Weisse  Max-Planck-Institut für Mathematik, Vivatsgasse 7, 53111 Bonn,
Germany  R. Gerstner  J. Sirker  Department of Physics and Astronomy and
Manitoba Quantum Institute, University of Manitoba, Winnipeg, Canada R3T 2N2

###### Abstract

We consider the spreading of a local operator $A$ in Euclidean time in one-
dimensional many-body systems with Hamiltonian $H$ by calculating the $k$-fold
commutator $[H,[H,[...,[H,A]]]]$. We derive general bounds for the operator
norm of this commutator in free and interacting fermionic systems with and
without disorder. We show, in particular, that in a localized system the norm
does grow at most exponentially and that the contributions of operators to the
total norm are exponentially suppressed with their length. We support our
general results by considering one specific example, the XXZ chain with random
magnetic fields. We solve the operator spreading in the XX case without
disorder exactly. For the Anderson and Aubry-André models we provide strict
upper bounds. We support our results by symbolic calculations of the
commutator up to high orders. For the XXX case with random magnetic fields,
these symbolic calculations show a growth of the operator norm faster than
exponential and consistent with the general bound for a non-localized system.
Also, there is no exponential decay of the contribution of operators as
function of their length. We conclude that there is no indication for a many-
body localization transition. Finally, we also discuss the differences between
the interacting and non-interacting cases when trying to perturbatively
transform the microscopic to an effective Hamiltonian of local conserved
charges by consecutive Schrieffer-Wolff transformations. We find that such an
approach is not well-defined in the interacting case because the
transformation generates $\sim 4^{\ell}$ terms connecting sites a distance
$\ell$ apart which can overwhelm the exponential decay with $\ell$ of the
amplitude of each individual term.

##  I Introduction

In a one-dimensional fermionic lattice model with short-range hoppings and
without interactions, any amount of disorder leads to a localization of the
eigenfunctions. This is known as Anderson localization [1, 2, 3, 4]. An open
question for the last 65 years has been what happens if interactions are
included. One of the difficulties in investigating this question is that
localization can no longer be defined in terms of single-particle
eigenfunctions; new indicators of localization are needed.

Studying the possibility of localization in the many-body case has gained
renewed interest in the last 17 years starting with the work by Basko,
Aleiner, and Altshuler [5]. This was followed by a large number of numerical
studies concentrating, in particular, on the Heisenberg (XXX) chain with
random magnetic fields [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].
By a Jordan-Wigner transform, this model is equivalent to a fermionic chain
with nearest-neighbor hoppings, density-density interactions, and random on-
site potentials. Based on an analysis of spectral properties such as the level
spacing as well as of dynamical properties such as the growth of the
entanglement entropy after a quantum quench from an unentangled initial state,
an apparent consensus was reached that this model shows a phase transition in
all eigenstates at a finite disorder strength, separating an ergodic from a
many-body localized phase. This consensus was first called into question by
studies showing that the transition point, defined for a finite system of size
$L$, appears to slowly shift to infinite disorder strength for $L\to\infty$
[20, 21]. These arguments were supported by studies of the number entropy
after a quantum quench [22, 23, 24, 25, 26, 27] and of the fidelity
susceptibility [28, 29].

Here we want to propose and follow a new avenue by studying the properties of
operators instead of wave functions. Our motivation is the following
observation: In an Anderson localized model, the projectors onto the
eigenstates $P_{n}=|E_{n}\rangle\langle E_{n}|$ are conserved quasi-local
operators, i.e., they are centered at a site in the lattice and decay
exponentially away from this site. These operators are also often called local
integrals of motion (LIOMs). Without loss of generality, we can order the
operators such that $P_{n}$ is an operator localized around lattice site $n$.
To be more precise, we can use the fermionic language in which the Anderson
Hamiltonian and any unitary transformation of this Hamiltonian is bilinear. We
can then write $P_{n}=\sum_{ij}a^{n}_{ij}c^{\dagger}_{i}c_{j}$ where
$c_{i}^{(\dagger)}$ is a fermionic annihilation (creation) operator at site
$i$ and $|a^{n}_{ij}|\sim\exp(|i-n|)\exp(|j-n|)$ are amplitudes. In
particular, the contribution of bilinear operators to $P_{n}$ is exponentially
suppressed with the distance they cover.

There is a unitary transformation $U$ which diagonalizes the microscopic, non-
interacting Anderson Hamiltonian $H$, bringing it into the form
$\tilde{H}=UHU^{-1}=\sum_{n}E_{n}P_{n}$ where $P_{n}$ is quasi-local and
conserved, $[H,P_{n}]=0$. Importantly, this unitary transformation itself is
quasi-local and will always map local microscopic operators $O_{i}^{\alpha}$
onto quasi-local operators, $\tilde{O}^{\alpha}_{i}=UO^{\alpha}_{i}U^{-1}$. As
we will show, this has important implications for the spreading of a local
operator in Euclidean time. We can write

$\displaystyle O_{i}^{\alpha}(\tau)$ $\displaystyle=$
$\displaystyle\text{e}^{\tau H}O^{\alpha}_{i}\text{e}^{-\tau H}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\underbrace{[H,[H,[...,[H,O_{i}^{\alpha}]]]]}_{\textrm{k-times}}\frac{\tau^{k}}{k!}.$

I.e., the spreading of the operator is encoded in the $k$-fold commutator
$[H,O_{i}^{\alpha}]^{(k)}$ with the Hamiltonian. It is important to note that
the Frobenius norm of this commutator is not affected by the unitary
transformation
$||[H,O_{i}^{\alpha}]^{(k)}||_{2}=||[\tilde{H},\tilde{O}_{i}^{\alpha}]^{(k)}||_{2}$.
For the latter we expect that the norm of operators in the commutator which
cover a distance of $l$ sites is exponentially suppressed with $l$ at fixed
order $k$. As we will show, this implies that the total norm of all terms in
the commutator can only grow exponentially with $k$.

In a many-body localized phase, it has been argued that—similar to the
Anderson case—a unitary transformation of the microscopic Hamiltonian $H$ to a
Hamiltonian $\tilde{H}$ of conserved quasi-local charges exists. The main
difference is assumed to be that these quasi-local charges $\tau^{z}_{n}$ are
interacting with each other [19, 17]

$\tilde{H}=\sum_{n}E_{n}\tau^{z}_{n}+\sum_{i,j}J_{ij}\tau^{z}_{i}\tau^{z}_{j}+\dots\,.$
(I.2)

Here, $J_{ij}$ are exponentially decaying coupling constants which describe a
slow dephasing between the conserved charges. For the Anderson case,
$J_{ij}=0$ and the $\tau^{z}_{n}$ can be chosen to be the projection operators
$P_{n}$. We note that the choice of the quasi-local conserved operators
$\tau^{z}_{n}$ is not unique. We will argue in the following that in the
interacting case a local operator in the microscopic model should also remain
local, up to exponential tails, under Euclidean time evolution as in the non-
interacting case. This implies, in particular, that the total norm in both
cases should at most grow exponentially with $k$.

This has to be contrasted with the general, non-localized case. In one
dimension, there is no phase transition in a clean system at finite
temperatures. Thus, the time evolution is always analytic and Eq. (I) implies
that the norm of the $k$-fold commutator $[H,O_{i}^{\alpha}]^{(k)}$ has to
grow slower than factorially. In Ref. [30], it has been shown that this is
indeed the case, however, the derived bound for the total norm grows much
faster than exponential and almost factorially. It has been hypothesized that
any generic quantum chaotic system will asymptotically saturate this bound
[30, 31]. Thus, we have a sharp contrast between the growth of the operator
norm of the $k$-fold commutator of a local operator with the Hamiltonian
between a localized system, where it can grow at most exponentially, and a
non-localized system, where an almost factorial growth is expected in general.
This qualitative difference in the operator norm growth is what we will
investigate here. To avoid misunderstandings, we want to stress already here
that what we mean by a local operator and what we will investigate in the
following is an operator $O^{\alpha}_{j}$ which resides on lattice site $j$ in
the microscopic model. In contrast, operators of the form
$O^{\alpha}=\sum_{j}O^{\alpha}_{j}$ are also often called local conserved
charges if $[H,O^{\alpha}]=0$ and are of particular importance in integrable
models. Operators of the type $O^{\alpha}$, which are sums of local densities
but which extend over the entire lattice, are also the type of operators whose
growth is investigated in Refs. [31, 32]. The latter type of operators are not
of interest for a study of localization and we will not consider them here.

In addition to analytical results, we will also present exact results from
symbolic calculations of commutators and norms. To be able to consider high
orders $k$ in the commutator, we have written code in Julia with Nemo [33]
which is able to handle symbolic computations with up to $\sim 200$ million
terms, limited only by the available memory. We note that we are able to
consider system sizes for nearest-neighbor models which are about twice as
large as the ones typically considered in exact diagonalizations.

Our paper is organized as follows: In Sec. II, we discuss general bounds for
the growth of the matrix norm $||[H,A]^{(k)}||$ for a local operator $A$ and a
nearest-neighbor Hamiltonian $H$. We review the result in the general
interacting case without localization [30] before discussing the free-fermion
case and the case of localization. In Sec. III, we then consider a specific
model, the $XXZ$ spin chain. We start with the free-fermion (XX) model without
disorder before discussing the Anderson case of quenched disorder. We also
consider the Aubry-André model which has a periodic potential incommensurate
with the lattice and shows a phase transition at a finite disorder strength
which, as we will demonstrate, is visible in the norm growth. We then turn to
the interacting Heisenberg (XXX) model and show that results obtained by
symbolic computations are inconsistent with localization for all studied
disorder strengths. In Sec. IV, we consider Schrieffer-Wolff transformations
which can be used to perturbatively transform the Hamiltonian into a basis of
quasi-local conserved operators. We show that such a perturbative construction
succeeds in the non-interacting case while it will typically fail once
interactions are included. The final section contains a brief summary of our
main results and a discussion of the remaining open issues.

##  II Bounds on operator norm growth

We want to consider how the norm of the $k$-th order commutator

$[H,A]^{(k)}\equiv\underbrace{[H,[H,[...,[H,A]]]]}_{\textrm{k-times}}$ (II.1)

between a one-dimensional lattice Hamiltonian

$H=\sum_{I}h_{I,I+1}\equiv\sum_{I}h_{I}$ (II.2)

which acts only locally at nearest-neighbor sites $I$ and $I+1$ and a one-site
operator $A$ grows with the order of the commutator $k$. We consider the
Hamiltonian as acting on $L$ lattice sites which therefore can be described by
an $M\times M$ matrix

$H=\sum_{I}\mathbbm{1}_{1}\otimes\cdots\otimes\mathbbm{1}_{I-1}\otimes
h_{I,I+1}\otimes\mathbbm{1}_{I+2}\otimes\cdots\otimes\mathbbm{1}_{L}$ (II.3)

with $M=n^{L}$ where $n$ is the number of degrees of freedom per site. We
place the operator $A$ in the middle of the chain and choose chain lengths
sufficiently large so that the support of $[H,A]^{(k)}$ is always less than
$L$.

We consider two entry-wise matrix norms for a matrix $A$ with elements
$a_{i,j}$:

$\displaystyle||A||_{1}$ $\displaystyle=$
$\displaystyle\sum_{i,j}|a_{i,j}|/\sqrt{\tr\mathbbm{1}_{M\times M}}$
$\displaystyle||A||_{2}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\tr(A^{\dagger}A)}{\tr\mathbbm{1}_{M\times
M}}}=\sqrt{\frac{\sum_{i,j}|a_{i,j}|^{2}}{\tr\mathbbm{1}_{M\times M}}}\,.$
(II.4)

We note that the $2$-norm is also called the Frobenius norm. For finite-
dimensional spaces it is also equivalent to the Hilbert-Schmidt norm which is
consistent with the standard scalar product in quantum mechanics at infinite
temperatures. We divide here by the constant factor $\tr\mathbbm{1}_{M\times
M}$ such that the identity operator acting on the entire system is normalized
to one. We note that the $2$\- norm is invariant under unitary
transformations. I.e., the operator growth of a local operator $A$ in the
microscopic model with Hamiltonian $H$ under Euclidean time evolution is
exactly the same as that of the transformed operator $\tilde{A}=UAU^{\dagger}$
in the effective model $H_{\textrm{eff}}=UHU^{\dagger}$ if $U$ describes a
unitary transformation. The $1$-norm, on the other hand, is not invariant
under basis changes. Let us therefore explain how we define the matrix $A$
corresponding to an operator. If $\\{\sigma_{i}^{\alpha}\\}$ is a
$f$-dimensional local operator basis, then we can always write

$A=\sum_{\\{\alpha\\}}a_{\alpha_{1},\dots,\alpha_{L}}\sigma^{\alpha_{1}}_{1}\otimes\cdots\otimes\sigma^{\alpha_{L}}_{L}\,.$
(II.5)

There are $f^{L}$ many coefficients $a_{\alpha_{1},\dots,\alpha_{L}}$, and the
$2$-norm is given by
$||A||_{2}=\sqrt{\sum_{\\{\alpha\\}}|a_{\alpha_{1},\dots,\alpha_{L}}|^{2}}$.
We can therefore consider as entries of our matrix $A$ the $f^{L}$ many
coefficients $\sqrt{\tr\mathbbm{1}_{M\times
M}}a_{\alpha_{1},\dots,\alpha_{L}}$. For the norms the order of these
coefficients does not matter. The $1$-norm is then given by
$||A||_{1}=\sum_{\\{\alpha\\}}|a_{\alpha_{1},\dots,\alpha_{L}}|$. In symbolic
calculations, where operators are obtained in the form (II.5), this means that
we simply have to sum up the absolute values, or the absolute values squared,
of the operator $A$ written in the local basis $\\{\sigma^{\alpha}_{i}\\}$.

The two matrix norms which we have defined in Eq. (II) are sub-additive and
the $1$-norm is an upper bound for the $2$-norm

$\displaystyle||A+B||_{p}$ $\displaystyle\leq$
$\displaystyle||A||_{p}+||B||_{p}$ $\displaystyle||A||_{2}$
$\displaystyle\leq$ $\displaystyle||A||_{1}\,.$ (II.6)

However, do to the rescaling by the factor $(\tr\mathbbm{1}_{M\times
M})^{-1/2}$ the two norms are no longer sub-multiplicative. We note, however,
that if
$A=J_{A}\sigma^{\alpha_{1}}_{1}\otimes\cdots\otimes\sigma^{\alpha_{L}}_{L}$
and $B=J_{B}\sigma^{\beta_{1}}_{1}\otimes\cdots\otimes\sigma^{\beta_{L}}_{L}$,
i.e., both operators have only one coefficient in Eq. (II.5) which is non-
zero, then

$||AB||_{p}=||A||_{p}||B||_{p}=|J_{A}||J_{B}|\,.$ (II.7)

While the $2$-norm (Hilbert-Schmidt norm) is physically more relevant because
it is induced by the standard scalar product in quantum mechanics, the
$1$-norm is often easier to calculate and provides an upper bound for the
$2$-norm.

###  II.1 General Case

The case of the Euclidean norm growth for a one-dimensional nearest-neighbor
model without any further restrictions on the form of the nearest-neighbor
Hamiltonian has been considered in Ref. [30]. Here we briefly recapitulate the
derivation of the general bound to contrast it with the tighter bounds which
we derive below in the special cases of free-fermion and localized models.

Let us define $H=\sum_{I}\sum_{a=1}^{C}h_{I}^{a}$ where $h_{I}^{a}$ are local
Hamiltonian densities of the form
$h_{I}^{a}=J^{a}_{I}\sigma_{I}^{a_{1}}\otimes\sigma_{I+1}^{a_{2}}$. We note
that the commutator (II.1) is zero if the local operator $A$ and the local
Hamiltonian densities $h^{a}_{I}$ do not form a connected cluster. Let
$\\{I_{1},\cdots,I_{k}\\}$ be a set of points which fulfills the adjacency
condition, then

$\displaystyle||[H,A]^{(k)}||$ $\displaystyle=$
$\displaystyle||\\!\\!\\!\sum_{\\{I_{1},\cdots,I_{k}\\}}\sum_{\\{a_{I}\\}}[h^{a_{I_{1}}}_{I_{1}},[h^{a_{I_{2}}}_{I_{2}},\cdots,[h^{a_{I_{k}}}_{I_{k}},A]]]||$
$\displaystyle\leq$
$\displaystyle\\!\\!\\!\sum_{\\{I_{1},\cdots,I_{k}\\}}\sum_{\\{a_{I}\\}}||[h^{a_{I_{1}}}_{I_{1}},[h^{a_{I_{2}}}_{I_{2}},\cdots,[h^{a_{I_{k}}}_{I_{k}},A]]]||$

where we have used the sub-additivity property (II). The commutator now is
determined by terms which are products of single operators in the
$\\{\sigma_{i}^{\alpha}\\}$ basis and we can use Eq. (II.7) implying that

$||h^{a_{I_{1}}}_{I_{1}}h^{a_{I_{2}}}_{I_{2}}\cdots
h^{a_{I_{k}}}_{I_{k}}A||=||h^{a_{I_{1}}}_{I_{1}}||\,||h^{a_{I_{2}}}_{I_{2}}||\cdots||h^{a_{I_{k}}}_{I_{k}}||\,||A||\,.$
(II.9)

There are $2^{k}$ terms of this kind. If we use again the sub-multiplicativity
property then we obtain

$\displaystyle||[H,A]^{(k)}||$ $\displaystyle\leq$
$\displaystyle\sum_{\\{I_{1},\cdots,I_{k}\\}}\sum_{\\{a_{I}\\}}\underbrace{||h||\cdots||h||}_{\text{$k$-many}}\,||A||$
(II.10) $\displaystyle\leq$
$\displaystyle(2CJ)^{k}||A||\sum_{\\{I_{1},\cdots,I_{k}\\}}$

where we have used that $\sum_{\\{a_{I}\\}}=C^{k}$ and
$||h^{a}_{I}||=J^{a}_{I}\leq J=\text{max}_{I,a}\\{J^{a}_{I}\\}$.

Eq. (II.10) is now a purely combinatorial problem. We need to figure out how
many different arrangements of the lattice points $\\{I_{1},\cdots,I_{k}\\}$
exist such that a connected cluster is formed. Without repeated indices, $j$
local adjacent Hamiltonians form a cluster with $j$ bonds. At each step, one
can either attach a local Hamiltonian to the right or to the left of the
existing cluster. Thus, as for the random walk in one dimension, there are
$2^{j}$ possibilities to form a cluster with $j$ bonds. Lastly, consider the
case that we have a cluster of length $j$ and $k\geq j$ local Hamiltonians
$h_{I}$. Then the question is how many ways there are to ‘distribute’ those
local Hamiltonians on the $j$ bonds of the cluster. This is the same question
as asking how many ways there are to partition $k$ elements into $j$ non-empty
sets. The answer is given by the Stirling numbers of the second kind,
$S(k,j)$. We therefore obtain

$\sum_{\\{I_{1},\cdots,I_{k}\\}}=\sum_{j=1}^{k}2^{j}S(k,j)=B_{k}(2)$ (II.11)

where $B_{k}(x)$ is the Bell (Touchard) polynomial. The general bound is
therefore given by

$s(k)\equiv||[H,A]^{(k)}||\leq||A||(2CJ)^{k}B_{k}(2)$ (II.12)

where $||\dots||$ can be either the $1$-norm or the $2$-norm defined earlier.
We note that this bound is also valid for any sub-multiplicative norm.

We are, in particular, interested in the asymptotic growth of the operator
norm for large $k$. In this limit, the Bell polynomial scales as

$B_{k}(2)\sim\frac{\exp[-k(1+\ln
W(k/2)-W^{-1}(k/2))]}{\sqrt{W(k/2)+1}}\text{e}^{-2}k^{k}$ (II.13)

where $W(x)$ is the Lambert W function [34]. Using the Stirling formula
$k^{k}\sim k!\exp(k)/\sqrt{2\pi k}$ we can rewrite the asymptotics as

$\displaystyle B_{k}(2)$ $\displaystyle\sim$
$\displaystyle\frac{\text{e}^{-2}k!}{\sqrt{2\pi k(W(k/2)+1})}$
$\displaystyle\times$ $\displaystyle\exp\left[k\frac{1-W(k/2)\ln
W(k/2)}{W(k/2)}\right].$

This means that the bound (II.12) grows faster than exponentially and almost
factorially. We note that the term in the exponential in the second line of
(II.1) only becomes negative for $k\gtrsim 20$ thus suppressing the initial
almost factorial growth. For large $k$, one can also use that $W(k)\sim\ln
k-\ln\ln k$ to see from Eq. (II.13) that $B_{k}(2)\sim(k/\ln k)^{k}$.

We can also ask the question how many distinct terms are possible in the
commutator $[H,A]^{(k)}$. Note that this question is different from the
question how many connected clusters exist which we addressed above because
connected clusters formed in different ways can result in the same operator.
Let us assume that our system has a basis of $f-1$ local operators plus the
local identity $\mathbb{I}_{j}$. At order $k$, clusters have at most length
$k+1$ and, if the operator $A$ was placed on site $j=0$, extend from sites
$-k,-k+1,\dots,0$ to the right. I.e., there are $k+1$ clusters of length $k+1$
with each shifted by one lattice site to the right compared to the previous
one. For the cluster starting at $-k$, we can put any of the $f-1$ local
operators or the identity on each of the $k+1$ sites leading to $f^{k+1}$
possible operators. For the cluster starting at site $k-1$, we have to put one
of the $f-1$ local operators at the right most site, $j=1$, to obtain
operators distinct from the ones already constructed. Thus, this cluster will
add $(f-1)f^{k}$ possible operators. The same construction also works for the
remaining $k-1$ clusters. We can therefore bound the total number of possible
distinct operators occuring in the commutator $[H,A]^{(k)}$ by

$n(k)\leq f^{k+1}+k(f-1)f^{k}=(kf-k+f)f^{k}\,.$ (II.15)

For a generic system, we thus expect that the number of distinct terms in the
commutator grows exponentially. If such a system does not have local operators
which are conserved, then the hypothesis is that the bound on the total norm
(II.12) is asymptotically tight [30]. We note again that by local we mean
operators acting on one or a few neighboring sites in the microscopic model,
not operators which are sums of local densities. Next, we consider non-
interacting systems where stricter, exponential bounds on the norm can be
derived.

###  II.2 Free Fermions

The main difference between a general, interacting system and a free-fermion
system is that in the latter case the commutator of a one-body operator with
the Hamiltonian always remains a one-body operator. Therefore many of the
connected clusters leading to the the bound (II.11) do not contribute. Here, a
free-fermion model is a system with a bilinear Hamiltonian of the form
$H=\sum_{j,k}h_{jk}c^{\dagger}_{j}c_{k}$ where $c_{j}^{(\dagger)}$ is a
fermionic annihilation (creation) operator at site $j$. We do not include pair
creation and annihilation terms in $H$ but note that the arguments in the
following remain valid qualitatively if they are included.

Our starting point to prove a bound in the free-fermion case is the last line
of Eq. (II.10). We now want to understand how many distinct connected clusters
there are which give a non-zero contribution to the sum. The Hamiltonian in
this case is, in general, a sum of $C$ one-body terms,
$H=\sum_{I}\sum_{a=1}^{C}h_{I}^{a}$ such as nearest-neighbor hoppings and
local potentials. Now consider a connected cluster
$\\{I\\}=\\{I_{1},\dots,I_{k}\\}$ which gives a non-zero contribution to the
sum in Eq. (II.10). Since the commutator $[H,A]^{(k)}$ is always a one-body
operator if $A$ is a one-body operator, there is one fermionic operator on a
site $i$ and another fermionic operator on a site $j$ for each such cluster.
Then, growing the cluster the only contributions which can potentially be non-
zero will come from adding a one-body term $h^{a}$ on one of the bonds
$I_{i-1},I_{i},I_{j},I_{j+1}$. I.e., there are at most $4$ possibilities to
build a new non-zero cluster with $k+1$ elements. While this is, in general,
overcounting the number of possibilities we can thus iteratively show that

$\sum_{\\{I_{1},\dots,I_{k}\\}}\leq 4^{k}$ (II.16)

resulting in the bound

$||[H,A]^{(k)}||\leq||A||(8CJ)^{k}\,.$ (II.17)

For a free-fermion system, the operator norm growth is thus at most
exponential [35].

We can also again consider the question how many distinct terms can at most
occur in the commutator $[H,A]^{(k)}$. As in the general case, the clusters
have length of at most $k+1$ and have $k+1$ different starting points
(leftmost lattice site in the cluster). The difference is that we now have to
distribute only two operators, a fermionic creation operator $c^{\dagger}_{i}$
and a fermionic annihilation operator $c_{j}$, over the cluster. For the first
cluster, we have $(k+1)$ possibilities to place the creation and $(k+1)$
possibilities to place the annihilation operator for a total of $(k+1)^{2}$
possible distinct terms. For the other clusters, we have to place either
$c^{\dagger}$ or $c$ on the rightmost site to get new distinct terms. The
other operator can then be placed on any of the $k+1$ sites for a total of
$2(k+1)-1=2k+1$ (placing both operators on the rightmost site results in only
one new term) distinct new terms for each of the $k$ remaining clusters. A
bound for the total number of distinct terms is thus given by

$n(k)\leq(k+1)^{2}+k(2k+1)=3k(k+1)+1\,.$ (II.18)

In a free-fermion system, the number of distinct operators in the commutator
$[H,A]^{(k)}$ thus grows quadratically with the order of the commutator $k$.

###  II.3 Localization

In a localized phase, a local operator $A$ is expected to remain localized to
a finite region of space of length $\xi_{\rm loc}$ (up to exponentially small
tails) when commuted $k$-times with the Hamiltonian. This means that we expect
that connected clusters of length $j\gg\xi_{\rm loc}$ give almost no
contribution to the norm. Neglecting such clusters leads to

$\sum_{\\{I_{1},\cdots,I_{k}\\}}\\!\\!\\!\\!\\!\\!\\!{\vphantom{\sum}}^{\prime}\;=\sum_{j=1}^{\min(k,\xi_{\rm
loc})}2^{j}S(k,j)\,.$ (II.19)

Asymptotically, for large $k$, we thus obtain

$\sum_{j=1}^{\xi_{\rm loc}}2^{j}S(k,j)\sim\sum_{j=1}^{\xi_{\rm
loc}}\frac{2^{j}j^{k}}{k!}<\frac{2^{\xi_{\rm loc}}}{(\xi_{\rm
loc}-1)!}\xi_{\rm loc}^{k}$ (II.20)

where we have used the asymptotic scaling $S(k,j)\sim j^{k}/j!$ of the
Stirling numbers of the second kind for $k\gg j$. The operator growth in a
localized phase is therefore at most exponential

$||[H,A]^{(k)}||\leq||A||\frac{2^{\xi_{\rm loc}}}{(\xi_{\rm
loc}-1)!}(2\xi_{\rm loc}CJ)^{k}\,.$ (II.21)

In deriving this bound, we have made the simplifying assumption that the
contributions of all clusters with lengths $j\gg\xi_{\rm loc}$ can be
neglected. Given that with increasing length also more possibilities exist to
distribute the local Hamiltonians $h^{a}_{I}$, one might worry that this
assumption is not justified. We will show in App. A that the argument remains
valid if clusters with $j>\xi_{\rm loc}$ are included.

We note that both the free-fermion case and the localized case show an
exponential growth of the operator norm. Thus, the total operator norm is
insufficient to distinguish a non-localized free-fermion phase from a
localized interacting or non-interacting phase. What does distinguish the two
phases though is how different terms in the commutator contribute to the norm.
If we denote by $[H,A]^{(k)}_{l}$ the terms in the commutator of order $k$
which have support on $l$ lattice sites, then we expect that in the localized
case

$s_{l}(k)\equiv||[H,A]^{(k)}_{l}||\sim\text{e}^{k}\text{e}^{-l}\,.$ (II.22)

For a non-localized phase, on the other hand, the initially local operator $A$
will spread over the entire lattice, implying that terms in the commutator
with $l\ll k$ should contribute approximately equally to the total norm.

To summarize, we hypothesize that the bound (II.12) is asymptotically tight
for an interacting, non-localized system which thus should have an operator
norm $s(k)$ which grows faster than exponential with the order of the
commutator $k$. In an interacting localized phase—the putative many-body
localized phase—the norm growth, as we have proven, can instead only be at
most exponential. For a free-fermion model, the norm growth is also always at
most exponential whether it is in a localized phase or not. In the latter
case, we can distinguish a localized from a non-localized phase by studying
the contributions $s_{l}(k)$ from terms which are supported on $l$ lattice
sites to the total norm of the commutator $s(k)$.

##  III The XXZ model with (quasi-)random fields

As a specific example, we want to consider the XXZ chain with random or quasi-
random magnetic fields

$H=\sum_{j}\left(\sigma^{x}_{j}\sigma^{x}_{j+1}+\sigma^{y}_{j}\sigma^{y}_{j+1}+\Delta\sigma^{z}_{j}\sigma^{z}_{j+1}+2h_{j}\sigma^{z}_{j}\right)\,.$
(III.1)

Here, $\sigma^{x,y,z}$ are Pauli spin-1/2 matrices, $\Delta$ characterizes the
spin exchange anisotropy, and $h_{j}$ are magnetic fields. As the operator
which we will commute with the Hamiltonian $H$ we consider a single local
$\sigma^{z}_{j}$. The model (III.1) is equivalent to a spinless fermion model
using the Jordan-Wigner transformation

$\displaystyle\sigma^{+}_{j}$ $\displaystyle=$ $\displaystyle
S_{j}c_{j}^{\dagger},\quad\sigma^{-}_{j}=S^{\dagger}_{j}c_{j},\quad\sigma^{z}_{j}=2c^{\dagger}_{j}c_{j}-1,$
$\displaystyle S_{j}$ $\displaystyle=$
$\displaystyle\exp\left[-i\pi\sum_{k=1}^{j-1}c_{k}^{\dagger}c_{k}\right]=\prod_{k=1}^{j-1}\left[1-2c_{k}^{\dagger}c_{k}\right]$
(III.2)

with $\sigma_{j}^{\pm}=(\sigma_{j}^{x}\pm i\sigma_{j}^{y})/2$ leading to

$\displaystyle\frac{H}{4}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{j}\left\\{c^{\dagger}_{j}c_{j+1}+c_{j+1}^{\dagger}c_{j}\right\\}$
$\displaystyle+$
$\displaystyle\sum_{j}\left\\{\Delta(n_{j}-1/2)(n_{j+1}-1/2)+h_{j}(n_{j}-1/2)\right\\}.$

Here $c_{j}^{(\dagger)}$ are fermionic annihilation (creation) operators and
$n_{j}=c_{j}^{\dagger}c_{j}$. In the following, we will investigate the non-
interacting case, $\Delta=0$, and the interacting, isotropic Heisenberg case,
$\Delta=1$.

###  III.1 Non-interacting case

We will first concentrate on the non-interacting free-fermion case,
$\Delta=0$, where localization is well understood. Here we want to show that
the general bounds for the total norm of the $k$-th order commutator derived
earlier do apply and that we can distinguish localized from non-localized
phases by resolving the contributions to the norm by the length of the spatial
support of each term in the commutator.

####  III.1.1 No disorder

The first specific case we will consider is the free-fermion case, $\Delta=0$,
without any magnetic fields, $h_{j}=0$. In this case, the entire structure of
the commutator $[H,\sigma^{z}_{0}]^{(k)}$ can be worked out analytically. We
find these results instructive and will discuss them in detail here to
motivate a description of the commutator by graphs where the vertices are the
terms generated and the edges denote their contribution to the 1-norm.

The full analytical solution is best understood in the fermionic language: the
commutator $[H,c_{0}^{\dagger}c_{0}]^{(k)}$ will only contain terms which are
one-particle operators, i.e., terms which contain exactly one annihilation and
one creation operator. More specifically, we find for the case that the order
of the commutator $k$ is odd that

$\displaystyle[H,\sigma_{0}^{z}]^{(k)}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\frac{k-1}{2}}\sum_{s=-\frac{k+1}{2}-l}^{\frac{k-1}{2}-l}a_{kls}\left\\{\sigma_{s}^{x}\bigotimes_{j=s+1}^{s+2l}\sigma_{j}^{z}\,\sigma_{s+2l+1}^{y}\right.$
$\displaystyle-$
$\displaystyle\left.\sigma_{s}^{y}\bigotimes_{j=s+1}^{s+2l}\sigma_{j}^{z}\,\sigma_{s+2l+1}^{x}\right\\}$
$\displaystyle=$
$\displaystyle-\frac{2}{i}\sum_{l=0}^{\frac{k-1}{2}}\sum_{s=-\frac{k+1}{2}-l}^{\frac{k-1}{2}-l}a_{kls}\left\\{c_{s}^{\dagger}c_{s+2l+1}-h.c.\right\\}$

with coefficients $a_{kls}$ which are discussed below. For $k$ even, we find
instead

$\displaystyle[H,\sigma_{0}^{z}]^{(k)}$ $\displaystyle=$
$\displaystyle\sum_{s=-\frac{k}{2}}^{\frac{k}{2}}a_{k0s}\,\sigma_{s}^{z}$
$\displaystyle+$
$\displaystyle\sum_{l=1}^{\frac{k}{2}}\sum_{s=-\frac{k}{2}-l}^{\frac{k}{2}-l}a_{kls}\left\\{\sigma_{s}^{x}\bigotimes_{j=s+1}^{s+2l-1}\sigma_{j}^{z}\,\sigma_{s+2l}^{x}\right.$
$\displaystyle+$
$\displaystyle\left.\sigma_{s}^{y}\bigotimes_{j=s+1}^{s+2l-1}\sigma_{j}^{z}\,\sigma_{s+2l}^{y}\right\\}$
$\displaystyle=$ $\displaystyle
2\sum_{s=-\frac{k}{2}}^{\frac{k}{2}}a_{k0s}\,(n_{s}-1/2)$ $\displaystyle-$
$\displaystyle
2\sum_{l=1}^{\frac{k}{2}}\sum_{s=-\frac{k}{2}-l}^{\frac{k}{2}-l}a_{kls}\left\\{c_{s}^{\dagger}c_{s+2l}+h.c.\right\\}.$

Note that the commutator has the form of a current density for $k$ odd and
that of a charge density for $k$ even. From the formulas (III.1.1) and
(III.1.1) we can immediately read off that there are
$n(k)=\frac{k+1}{2}(k+1)\cdot 2=k+1+\frac{k}{2}(k+1)\cdot 2=(k+1)^{2}$
different terms generated at order $k$, both in the odd and in the even case.
We confirm this result by calculating the commutator using symbolic
manipulations up to high order, see Fig. 1.

Figure 1: XX model: Without disorder, the number of distinct terms $n(k)$ in
$[H,\sigma^{z}_{0}]^{(k)}$ grows as $(k+1)^{2}$ while it grows according to
Eq. (III.13) with disorder. The general bound (II.18), valid for any free-
fermion model, is also shown. Figure 2: In the XX model without disorder the
norms grow as $s^{(1)}(k)=8^{k}$ and $s^{(2)}(k)=2^{k}\binom{2k}{k}$. Note
that on this scale both results are very close to each other. The $2$-norm is
only slightly smaller than the $1$-norm. Also shown is the general bound, Eq.
(II.17), with $C=2$.

We want to stress that the quadratic growth of the number of terms is a
special property of a free fermionic model where a one-body operator, when
commuted with the Hamiltonian, always remains a one-body operator. What
changes is only the separation between the creation and annihilation operator
(support of the one-body operator) and the number of different sites where the
first of the fermionic operators is located. Both the allowed separations as
well as the possible sites for the first fermionic operator grow linearly in
$k$, leading to the overall quadratic growth in the number of distinct terms
in the commutator. We already want to note here that this behavior is expected
to change dramatically in the interacting case.

In the simple, non-disordered XX case considered here, we can even take it one
step further and calculate the coefficients $a_{kls}$ exactly. These
coefficients are completely determined by the coefficient matrices $C_{k}$
which are defined by

$(C_{k})_{ij}=2^{k}{{k}\choose{i-1}}{{k}\choose{j-1}}$ (III.6)

with $1\leq i,j\leq k+1$. In the following, we concentrate on the case of $k$
odd which has a slightly simpler structure, see Eq. (III.1.1). In this case we
find, in particular, that

$a_{kls}=\text{i}(-1)^{s}(C_{k})_{s+\frac{k+3}{2}+l,\frac{k+1}{2}-l}.$ (III.7)

From the result (III.7) we can now determine the total norms of the commutator
and, furthermore, we can even determine which operators contribute how much to
the total norm. From Eq. (III.1.1) and Eq. (III.7) we find that the set of the
coefficients $\\{|a_{kls}|\\}$ is the same as $\\{(C_{k})_{ij}\\}^{2}$, i.e.,
each entry in $C_{k}$ appears two times. We therefore find

$s^{(1)}(k)\equiv||[H,\sigma^{z}_{0}]^{(k)}||_{1}=2\sum_{i=1}^{k+1}\sum_{j=1}^{\frac{k+1}{2}}(C_{k})_{ij}=8^{k}$
(III.8)

where we have used the explicit form of the coefficient matrix (III.6).
Similarly, we find

$s^{(2)}(k)\equiv||[H,\sigma^{z}_{0}]^{(k)}||^{2}_{2}=2\sum_{i=1}^{k+1}\sum_{j=1}^{\frac{k+1}{2}}(C_{k})_{ij}^{2}=2^{2k}{{2k}\choose{k}}^{2}.$
(III.9)

Asymptotically, the $2$-norm therefore grows as

$s^{(2)}(k)=2^{k}{{2k}\choose{k}}\approx\frac{8^{k}}{\sqrt{\pi k}}$ (III.10)

where we have used Stirling’s formula. This shows that the $2$-norm grows also
exponentially in $k$ but with a correction which makes it always smaller than
the $1$-norm as expected on general grounds, see Eq. (II). We can again check
these results for finite $k$ by symbolic manipulations, see Fig. 2.

Finally, we can study how much terms which extend over a distance $l$
contribute to the total norm. We consider the case $k$ odd where the amplitude
$a_{kls}$ for terms which have support on $2l+2$ sites is given by Eq.
(III.7). We find again that the same entries appear multiple times and that
all prefactors for a given $l$ are represented by the column
$j=\frac{k+1}{2}-l$ of the coefficient matrix $(C_{k})_{ij}$. Summing over all
initial sites we thus obtain

$\displaystyle s^{(1)}_{2l+2}(k)$ $\displaystyle=$
$\displaystyle||[H,\sigma^{z}_{0}]^{(k)}_{2l+2}||_{1}=\sum_{i}|a_{kli}|$
$\displaystyle=$ $\displaystyle
2\sum_{i=1}^{k+1}(C_{k})_{i,\frac{k+1}{2}-l}=2^{2k+1}{{k}\choose{\frac{k-1}{2}-l}}$
$\displaystyle\rightarrow$ $\displaystyle
s^{(1)}_{l}(k)=2^{2k+1}{{k}\choose{\frac{k+1-l}{2}}}$

where in the last line we transformed the equation to show the support on $l$
sites rather than on $2l+2$ sites. In the limit $k\gg l$ we can apply
Stirling’s formula and find

$s^{(1)}_{l}(k)\approx\sqrt{\frac{8}{k\pi}}8^{k}\,.$ (III.12)

This means that all the terms contribute equally independent of the length of
their support $l$ which clearly indicates delocalization. This analytical
result can be checked by symbolic manipulations, see Fig. 3.

Figure 3: For $k\gg l$ the norm of operators with support on $l$ sites is
asymptotically described by Eq. (III.12) and is independent of $l$.

While we have been able to exactly calculate the $k$-th order commutator in
the non-interacting, non-disordered case, see Eqs. (III.1.1, III.1.1, III.7),
this is clearly no longer possible, even in principle, if random fields are
included. On a positive note, this is also more than what we need. To
distinguish between free-fermionic localized and non-localized phases and
interacting localized and non-localized phases it is sufficient to consider
the total commutator norm and the norm of the terms in the commutator which
are supported on $l$ sites. This can be achieved more efficiently by
developing graphs where the vertices represent the possible terms at order $k$
and the edges represent the contribution to the $1$-norm. Fig. 4 shows such a
graph for the non-interacting model. The case without random magnetic fields
discussed here corresponds to the case $D=0$ where only the terms in the left
column are present.

$\sigma^{z}$$\sigma^{x}\sigma^{y}$$\sigma^{x}\sigma^{z}\sigma^{x}$$\sigma^{x}\sigma^{z}\sigma^{z}\sigma^{y}$$\sigma^{x}\sigma^{z}\sigma^{z}\sigma^{z}\sigma^{x}$$\sigma^{x}\sigma^{x}$$\sigma^{x}\sigma^{z}\sigma^{y}$$\sigma^{x}\sigma^{z}\sigma^{z}\sigma^{x}$$\sigma^{x}\sigma^{z}\sigma^{z}\sigma^{z}\sigma^{y}$$8$$4$$4$$4$$4$$4$$4$$4$$8D$$8D$$8D$$8D$

⋮ ⋮

Figure 4: Graph for the commutator $[H,\sigma^{z}_{0}]^{(k)}$ in the non-
interacting case, $\Delta=0$. The vertices represent the possible terms with
support on $l$ sites ($l=1,2,\dots$ from top to bottom) while the edges
represent the contribution to the $1$-norm. The number of traversed edges
equals the order of the commutator $k$. Terms obtained by
$\sigma^{x}\leftrightarrow\sigma^{y}$ are grouped together into a single
vertex. Without disorder, $D=0$, only the terms in the left column will be
present. Figure 5: $s^{(1)}(k)$ and $s^{(2)}(k)$ for the Anderson model with
$D=2$ and $D=10$. The results are averaged over $100$ samples. Note that the
$2$-norm is only slightly smaller than the $1$-norm and the two norms are
almost indistinguishable on this scale. The bound (III.14) is shown as well.

To obtain $s^{(1)}(k)$, we have to sum over all weighted walks through the
graph which start at the top and have length $k$ by multiplying the values of
all traversed edges, irrespective of the vertex the walk ends on. If, on the
other hand, we want to calculate only the contributions of terms which have
support on $l$ sites, $s^{(1)}_{l}(k)$, then we have to sum only over those
weighted walks with $k$ steps which end at a vertex with terms on $l$ sites.
We discuss this approach to calculating or bounding the norm of the commutator
in more detail in App. B. Fig. 4 also already shows how the graph has to be
modified in the Anderson case which we discuss next.

####  III.1.2 Anderson Localization

Next, we consider the case of quenched disorder where the magnetic fields
$h_{j}$ in the Hamiltonian (III.1) are drawn randomly from a box distribution,
$h_{j}\in[-D,D]$. We note that in the one-dimensional case considered here any
amount of disorder is a relevant perturbation. The model is localized for all
$D>0$. Since $h_{j}$ is random, we cannot expect exact results for the norm.
Instead, we will derive strict upper bounds.

What we can do exactly, however, is to count the number of distinct terms
$n(k)$ in the commutator which occur at order $k$ and which is independent of
the specific disorder realization. This is easiest done by making use of the
graph Fig. 4 and by taking also into account the different possibilities where
in the lattice these terms can be located. Doing so we find

$n(k)=2k(k+1)+\left\\{\begin{array}[]{cc}1&\mbox{even}\\\
1-k&\mbox{odd.}\end{array}\right.$ (III.13)

This result is compared to symbolic calculations in Fig. 1. We stress once
more that Eq. (III.13) is independent of the type of disorder or the specific
disorder configuration and is thus also valid in the Aubry-André case
discussed in Sec. III.1.3. Note that while the precise number of terms has
changed as compared to the non-disordered case, importantly the number of
terms still only increases quadratically with $k$.

Figure 6: Left column: Norm ratio $s^{(1)}_{l}(k)/s^{(1)}_{1}(k)$ for a fixed
$k$ as a function of $l$ for the Anderson model with $D=2$ (top row) and
$D=10$ (bottom row). Right column: $s^{(1)}_{l}(k)$ for a fixed $l$ as a
function of $k$. The data are averaged over $100$ samples.

To bound the $1$-norm, we use the graph shown in Fig. 4. With random fields,
the additional terms in the right column and the corresponding additional
paths are present. In the first step, we collect a factor of $8$ in the
$1$-norm. In every other step, we can collect at most a factor $8+8D$ where we
have bounded the disorder by $|h_{j}|\leq D$. We therefore obtain the
following bound

$s^{(1)}(k)\leq 8(8+8D)^{k-1}=8^{k}(D+1)^{k-1}.$ (III.14)

As expected, we find that the norm has an exponential bound. Note that this
bound increases with increasing disorder $D$, see Fig. 5. Also note that for
$D=0$ we obtain the exact result for the norm in the non-disordered case, Eq.
(III.8).

To see the difference between a free-fermionic model, where an initially
localized operator spreads over the entire chain, and a system where the
operator remains localized it is thus not sufficient to study the total
operator norm. Instead, we have to investigate the spatial structure of the
terms in the commutator at order $k$. Based on the graph in Fig. 4, we can
derive a bound for the contribution of terms with support on $l$ sites to the
norm at order $k$. We use that $|h_{j}|\leq D$. The derivation of this bound
is discussed in detail in App. B.2. Here, we will consider the simplest
possible approach where we only consider those terms with support on $l$ sites
which have the largest possible number of factors $h_{j}$. This approach is
justified for $D\gg 1$ and, as we will argue, gives the same asymptotics in
the limit $k\gg l\gg 1$ as the more sophisticated approach discussed in App.
B.2. Let us first consider terms which have support on a single site, i.e.,
which are of the form $\sigma^{z}_{j}$. As is obvious from the graph in Fig.
4, such terms are only generated at even orders of $k$. We collect the largest
number of factors $h_{j}$ if we go down one level, then travel $k-2$ times
along the vertical direction, and then go back up to the top of the graph. The
$1$-norm of such terms is thus bounded by

$s^{(1)}_{1}(k)\leq 8\times(8D)^{k-2}\times 4=2^{2k+1}(2D)^{k-2}\,.$ (III.15)

Next, we consider the norm of terms with support on $l\geq 2$ sites with the
largest number of possible $h_{j}$ factors. In the first step, we have to go
one level down in the graph in Fig. 4 which contributes a factor of $8$. To
get to terms with support on $l\geq 2$ sites we have to take another $l-2$
steps down, contributing a factor $4^{l-2}$. We can distribute these steps
among a total of $k-1$ remaining steps, i.e., there are $\binom{k-1}{l-2}$
possibilities. Finally, we get a factor $(8D)^{k-1-(l-2)}$ from the remaining
steps along the vertical lines in the graph. Thus, the one norm of such terms
is bounded by

$s^{(1)}_{l}(k)\leq\binom{k-1}{l-2}2^{2k+1}(2D)^{k-l+1}\,.$ (III.16)

Figure 7: Anderson model with $D=10$. Left: Norm ratio
$s^{(1)}_{l}(k)/s^{(1)}_{3}(k)$ for a fixed $l$ as a function of $k$. Right:
The norm ratio $s^{(1)}_{l}(k\to\infty)/s^{(1)}_{3}(k\to\infty)$ decays
exponentially with the support $l$. The data are averaged over $100$ samples.

The relative contribution to the norm of terms with support on $l$ sites when
normalized to the contribution coming from single-site operators is thus given
by

$\frac{s_{l}(k)}{s_{1}(k)}\lesssim\binom{k-1}{l-2}(2D)^{3-l}\approx\frac{1}{\sqrt{2\pi
l}}\left(\frac{ke}{l}\right)^{l}e^{-\frac{l^{2}}{2k}}(2D)^{3-l}\,.$ (III.17)

We note that while for a fixed $k$ we indeed find that the ratio (III.17) goes
exponentially to zero for large $l$, there is a maximum in the ratio which
shifts to larger $l$ with increasing $k$. This is not quite the behavior we
expect. Instead, we expect that with increasing $k$ the ratio (III.17)
converges to a function which is exponentially decreasing in $l$. Here it is
important to note that (III.17) is only an upper bound. In particular, we have
only considered those paths which have a maximum number of $h_{j}$ factors.
There are, however, many additional paths which can also interefere with each
other and we expect that the prefactor in (III.17) will eventually converge
with $k$ for a given $l$. This is consistent with what we observe from
symbolic manipulations to finite orders of $k$, see Fig. 6. We note that in
contrast to the non-disordered case where $s_{l}(k)$ grows exponentially and
independent of $l$ for $k\gg l$, see Eq. (III.12) and Fig. 3, we find that
with quenched disorder

$s_{l}(k)\sim\text{e}^{k}D^{-l}\,,$ (III.18)

see Fig. 7. The spatial structure of the norm contributions $s_{l}(k)$ thus
provide a clear distinction between the non-localized and localized phases.

####  III.1.3 Aubry-André model

To further check that the exponential decay of
$s_{l}(k)=||[H,\sigma^{z}_{0}]^{(k)}_{l}||$ with $l$ is indeed the signature
of a localized system, we consider the Aubry-André model next. This is a model
which is typically formulated as a fermionic model (III) where instead of
quenched disorder we have a periodic potential

$h_{j}=\frac{\lambda}{2}\cos(2\pi\beta j+\phi)\,.$ (III.19)

Then, for $\beta$ Diophantine and almost any $\phi$—excluding a set of measure
zero of resonant phases—the eigenstates of the Aubry-André model are localized
for $\lambda>2$ [36, 37]. A common choice is to use the golden ratio
$\beta=(1+\sqrt{5})/2$. The model thus provides an important test case where
both a non-localized and a localized phase exist. Below the transition in the
non-localized phase we expect that $s_{l}(k)$ is qualitatively described by
Eq. (III.12), i.e., grows with $k$ independent of $l$ for $k\gg l$. In the
localized phase, on the other hand, we expect that $s_{l}(k)$ is exponentially
localized, see Eq. (III.18). Note that the graph for the Aubry-André model is
the same as the one in the Anderson case shown in Fig. 4 with $D\to\lambda/2$
on the edges. With this replacement, also the bound (III.14) holds. We note
that the considerations which let us to Eq. (III.18) in the Anderson case were
valid only in the case of strong disorder. In this limit, they remain valid
for the Aubry-André model. That the exponential decaying structure for the
spatial support of the norm remains unchanged for all $D\neq 0$ in the
Anderson case while there is a transition to the structure (III.12), which is
independent of $l$ for $k\gg l$, for the Aubry-André model is beyond the
strong disorder arguments we used in the previous section. All these results
and expectations outlined above are fully consistent with the symbolic
computations shown in Fig. 8.

Figure 8: Aubry-André model for $\beta=(1+\sqrt{5})/2$, $\phi=0$ where
$\lambda=0.5$ (top row) and $\lambda=10$ (bottom row). For $\lambda=0.5$,
$s^{(1)}_{l}(k)$ has the same exponential asymptotics independent of $l$
clearly showing that the model is not in the localized phase. In contrast, for
$\lambda=10$ the norm scales as $s^{(1)}_{l}(k)\sim\text{e}^{k}\text{e}^{-l}$,
i.e., is exponentially suppressed with increasing $l$ for $k$ fixed,
demonstrating that the model is now in the localized phase.

While other methods might be better suited to determine the exact phase
transition point, there is a clear qualitative difference in the scaling of
$s_{l}(k)$ between the extended and the localized phase. The operator growth
$s_{l}(k)$ therefore allows for an unambiguous distinction between these two
phases.

###  III.2 Interacting case

With the results we have established in the non-interacting non-localized and
localized cases, we can now look at the interacting isotropic Heisenberg (XXX)
model, i.e., the Hamiltonian (III.1) with $\Delta=1$. A fundamental general
difference in interacting as compared to non-interacting models is that one-
body operators do not remain one-body operators when commuted with the
Hamiltonian. This leads to a qualitatively different phenomenology in the
operator growth with exponentially many more terms and paths.

We also note that the Heisenberg model without disorder is a so-called
integrable model which has, what is commonly called, an infinite set of local
conserved charges $O$ with $[H,O]=0$ in the thermodynamic limit. It is crucial
to note though that ‘local’ is meant here in an entirely different way. These
charges are in no way localized but rather can be written as $O=\sum_{j}o_{j}$
where $o_{j}$ is a local density acting on a finite number of lattice sites.
I.e., the charge densities are local but the entire conserved operator $O$ is
not and rather acts on the entire lattice. The trivial examples are the
Hamiltonian itself, $H=\sum_{I}h_{I}$ and the magnetization
$M=\sum_{I}\sigma^{z}_{I}$. What is special about integrable models is that an
infinite set of operators of this type exist which allows to uniquely
characterize all eigenstates by the corresponding quantum numbers. The
eigenspectrum of the model can thus be calculated exactly by using the Bethe
ansatz [38, 39]. For the question of local operator growth, integrability is
not relevant and we do not expect to see any difference to a non-integrable
model. Even if we break integrability, the magnetization $M$ is still
conserved as long as $U(1)$ symmetry is preserved. I.e., we are in any case
studying the operator growth of a local operator $\sigma^{z}_{j}$ which
belongs to a globally conserved charge. Differences between integrable and
non-integrable model are, of course, expected if we study the Euclidean time
dynamics of operators of the type $O=\sum_{j}o_{j}$ [31]. However, for the
Euclidean time dynamics of local operators in a generic system, the only thing
which matters is whether or not this system has local conserved charges, i.e.,
whether it is localized or not. We also note that integrability is immediately
destroyed if we introduce even an infinitesimal amount of disorder.

We will first briefly consider the case without disorder before turning to the
most interesting case, the XXX model with random magnetic fields
$h_{j}\in[-D,D]$ which has been suggested to contain a localized phase for
sufficiently strong disorder.

####  III.2.1 No disorder

From the exact solution of the model, we know for sure that the eigenstates of
the XXX model are not localized. We can calculate the number of distinct terms
$n(k)$ in the commutator $[H,\sigma_{0}^{z}]^{(k)}$ at order $k$ by symbolic
calculations up to order $k=18$ exactly. We have not been able to identify the
series but asymptotically the number of terms $n(k)$ does grow exponentially
with a multiplicative algebraic correction which is consistent with the
general bound (II.15) which we derived earlier. The data and a corresponding
fit are shown in Fig. 9.

Figure 9: In the interacting model with $\Delta=1$, the number of different
terms $n(k)$ in the commutator $[H,\sigma_{0}^{z}]^{(k)}$ grows exponentially
(symbols) and is very well fitted by $n(k)\sim 1.514\,k^{0.803}\exp(0.922k)$
in the non-disordered case and $n(k)\sim 1.492\,k^{0.815}\exp(0.998k)$ in the
case with disorder (solid lines). The general bound (II.15) with $f=3$ and
$f=4$ is shown as well.

The fit works extremely well with a least square deviation of $\chi^{2}\sim
10^{-2}$. Also shown in Fig. 9 is the general bound on $n(k)$, Eq. (II.15).
The local operator basis for the Heisenberg model has dimension $f=4$ and
consists of the operators ${\mathbbm{1},\sigma^{x},\sigma^{y},\sigma^{z}}$. We
note, however, that only a small subset of terms in $[H,\sigma^{z}_{0}]^{(k)}$
contains an identity in between Pauli matrices. Ignoring such terms, we can
set $f=3$ in Eq. (II.15) which gives a much tighter bound, see Fig. 9.

We note that the number of distinct terms is different and always smaller than
the number of connected clusters which we have calculated in Eq. (II.11). The
reason is that different connected clusters can result in the same term. For
the calculation of the norm this does not matter though. We can bound the
contribution of each connected cluster as in Eq. (II.12) even if some of the
resulting terms are the same. In Fig. 10 we show that this general bound
indeed describes the norm growth well if we use the bound on the local
Hamiltonian, $||h^{a}_{I}||\leq J$, as a fitting parameter.

Figure 10: The 1-norm in the interacting case with $\Delta=1$ grows faster
than exponential for all disorder strengths shown, consistent with the
rescaled general bound, Eq. (II.12). The data are averaged over $30$ disorder
realizations for $D=2.5$, $D=5$ and $D=12.5$ and over 50 realizations for
$D=25$.

In particular, the norm clearly grows faster than exponential for the orders
of the commutator which we can handle symbolically. We also note that we
obtain qualitatively similar results if we use the $2$-norm instead of the
$1$-norm, see Fig. 11.

Figure 11: The $2$-norm ratio $s^{(2)}(k+1)/s^{(2)}(k)$ for $D=0,5,12.5,25$.
The data for $D=5$ and $D=12.5$ are averaged over $30$, and the data for
$D=25$ over $50$ realizations. The norm ratio has to become a constant for
large $k$ if the norm grows exponentially.

This difference between the non-interacting and the interacting case becomes
even more evident if we consider the norm $s^{(1)}_{l}(k)$ of terms supported
on $l$ sites, see Fig. 12.

Figure 12: The norm of operators $s^{(1)}_{l}(k)/s^{(1)}_{1}(k)$ with support
on $l$ sites for the Heisenberg model with disorder $D=0,5,12.5,25$. The data
for $D=5$ and $D=12.5$ are averaged over $30$, the data for $D=25$ over $50$
disorder realizations. There is no indication for a convergence with
increasing order $k$ to an exponential decay with $l$ even for very strong
disorder. For $D=0$ and $D=25$ the results for the non-interacting case are
shown for comparison (symbols with dotted lines).

Whereas in the non-interacting, non-disordered case the norm of the terms
supported on $l$ sites grows asymptotically in the same way independent of $l$
if the order of the commutator $k$ is much larger than $l$, see Eq. (III.12),
$s^{(1)}_{l}(k)/s^{(1)}_{1}(k)$ has a maximum which shifts to larger and
larger $l$ with increasing $k$ in the interacting case. For large $k$, the
total norm is thus completely dominated by terms which have support on a large
number of lattice sites. This is to be expected, given that the number of
distinct terms in the commutator is increasing exponentially which also leads
to a rapidly growing number of possible paths between those terms, see also
the graph in Fig. 13.

$\sigma^{z}$$\sigma^{x}\sigma^{y}$$\sigma^{x}\sigma^{z}\sigma^{x}$$\sigma^{x}\sigma^{x}$$\sigma^{x}\sigma^{z}\sigma^{y}$$\sigma^{z}\sigma^{x}\sigma^{x}$$\sigma^{z}\sigma^{x}\sigma^{y}$$\sigma^{x}\sigma^{y}\sigma^{x}\sigma^{x}$$\sigma^{x}\sigma^{x}\sigma^{x}\sigma^{y}$$\sigma^{z}\sigma^{x}\sigma^{y}\sigma^{x}$$\sigma^{z}\sigma^{x}\sigma^{z}\sigma^{y}$$\sigma^{x}I\sigma^{y}$$\sigma^{x}\sigma^{z}\sigma^{z}\sigma^{y}$$\sigma^{z}\sigma^{x}\sigma^{z}\sigma^{x}$$\sigma^{x}\sigma^{z}\sigma^{z}\sigma^{x}$$\sigma^{x}\sigma^{y}\sigma^{y}\sigma^{x}$$\sigma^{x}\sigma^{y}\sigma^{x}\sigma^{y}$$\sigma^{x}I\sigma^{x}$$\sigma^{z}\sigma^{z}$$\sigma^{z}I\sigma^{z}$$8$$4$$4$$4$$8D$$8D$$4$$2$$4$$8D$$2$$2$$2$$2$$2$$4$$4$$4$$4$$4$$4$$2$$2$$2$$2$$2$$2$$2$$2$
Figure 13: Graph for the $1$-norm $s^{(1)}(k)$ up to order $k=3$ for the
Heisenberg model with random fields $h_{j}\in[-D,D]$. The number of terms
grows exponentially with the length of their support.

####  III.2.2 Many-body localization

From the perspective of operator spreading, localization in a many-body system
is thus an extremely challenging requirement. In a graph with an equal number
of edges at each level as realized in a non-interacting model, see Fig. 4, it
is easy to see—at least in the limit of large disorder—how the norm of terms
with support on only a few sites can dominate the total norm. The exponential
growth of edges in the interacting case (see Fig. 13), on the other hand,
suggests that the total norm for large $k$ will be entirely dominated by terms
with support on many lattice sites. I.e., a local operator is generically
expected to delocalize. This aspect also makes it difficult to see how a
perturbative construction of local integrals of motion can be arranged. This
is a point which seems to have not been sufficiently appreciated in Ref. [40,
41] where such a construction was suggested. We will get back to this point in
Sec. IV. The only way a local operator can avoid to spread through the entire
system in the many-body case thus seems to be intricate interference effects
between the exponentially many different paths in the graph Fig. 13. We cannot
prove that such interference effects do not exist. Instead, we will consider
data from symbolic computations up to order $k=16$ and show that there is no
evidence for localization. We note that at order $k$, we are considering
operators which are spread over $2k+1$ sites in the lattice. I.e., at order
$k=16$ we are dealing effectively with a chain of $33$ lattice sites which is
well beyond what is possible in exact diagonalizations. We will also consider
disorder strengths up to $D=25$ which is much larger than the critical values
which have been suggested in the past [8, 10, 11, 42].

Let us start by considering the number of distinct terms in the commutator at
order $k$ which is shown in Fig. 9. We find that the number of terms grows
faster in the disordered than in the case without disorder. This is easy to
understand: additional terms appear in the commutator due to the random
magnetic field terms while all the terms present before remain. The different
classes of terms which appear up to order $k=3$ are shown in the graph in Fig.
13. Next, we consider the growth of the total norm which should be at most
exponential in a localized phase as we have proven in Eq. (II.21), see also
App. A. The data, however, are overall much better described by the general
bound (II.12) than by an exponential, with the norm bound of the local
Hamiltonian, $||h^{a}_{I}||\leq J$, used as a fitting parameter. While this
does not exclude that for larger $k$ the scaling eventually does become
exponential, we can say with confidence that no indications for localization
are observed in spin chains up to $33$ sites. These system sizes are much
larger than those used in exact diagonalizations, the method which claims of a
many-body localized phase have been mostly based on [8, 10, 11, 19]. The
stronger than exponential growth of the norm is further supported by
considering the ratio of norms $s^{(2)}(k+1)/s^{(2)}(k)$ which is shown in
Fig. 11. For the orders of the commutator which we can handle symbolically,
there is no indication that this norm ratio saturates, even for extremely
strong disorder, as it has to in a localized case where $s(k)\sim\exp(k)$. As
the figure shows, the behavior in the interacting case is in stark contrast to
the localized non-interacting case.

Finally, let us also consider the norm ratio $s^{(1)}_{l}(k)/s^{(1)}_{1}(k)$
of terms with support on $l$ sites which is shown in Fig. 12 for various
disorder strengths. Even for very strong disorder $D=25$ there is no
indication of a convergence to an exponential decay with $l$. To the contrary,
there is always a maximum which shifts to larger $l$ with increasing order $k$
just as in the case without disorder.

Figure 14: Norm ratio $s^{(1)}_{l}(k)/s^{(1)}_{3}(k)$ for a fixed $l$ as a
function of $k$ for the Heisenberg model with $D=12.5$ (top) and $D=50$
(bottom). These results are to be compared with the non-interacting case shown
in Fig. 7.

To further highlight this point, we show in Fig. 14 the norm ratio
$s^{(1)}_{l}(k)/s^{(1)}_{3}(k)$ as a function of the order $k$ for various
fixed $l$. In contrast to the non-interacting case shown in Fig. 7, the curves
do not saturate for large $k$ but rather are consistent with the weight in the
total norm shifting to larger and larger $l$ with increasing $k$. While we
cannot exclude that this trend changes for even larger $k$ than the ones we
can handle in symbolic calculations, we want to stress once more that the
terms in the commutator are spread over $33$ sites for $k=16$. This is already
much larger than the system sizes typically considered in exact
diagonalizations. At the very least, we can say that for chains of lengths up
to $33$ sites, there is no indication that a many-body localized phase will be
realized even at disorder strengths up to $D=25$.

##  IV Schrieffer-Wolff transformations

Figure 15: Left: The number of terms in the Hamiltonian with support on $l$
sites after performing the Schrieffer-Wolff transformation (IV.9) up to order
$k$ for the XXX chain with $D=5$ and $N=20$ lattice sites. For large $k$, this
number increases approximately exponentially with $l$. Right: The
corresponding sum of $1$-norms of terms with support on $l$ sites averaged
over $30$ realizations. While each term is expected to have an amplitude $\sim
D^{-2l}$ if there are no resonances, the exponentially increasing number of
terms means that the total coupling of distant sites is no longer controlled
by a small parameter. For example, we have $s_{1}^{(1)}(k)\approx
s_{7}^{(1)}(k)$ for $k=4$.

We have so far considered the commutator of a microscopic lattice Hamiltonian
$H$ with a local operator $\sigma^{z}_{i}$. If this Hamiltonian is in a
localized phase, then it is unitarily equivalent to a Hamiltonian of
(interacting) conserved charges $\tau^{z}_{j}$ which are localized around
lattice site $j$, see Eq. (I.2). Our argument has been that the local operator
$\sigma^{z}_{i}$ then only has significant overlap with conserved charges
$\tau^{z}_{j}$ which are localized near site $i$. Thus, the operator
$\sigma^{z}_{i}$ will remain localized—up to exponential tails—under Euclidean
time evolution.

An alternative but related approach is to try to explicitly construct the
unitarily equivalent Hamiltonian written in terms of local conserved charges
$\tau^{z}_{j}$ [40, 41, 43, 44, 45, 29]. It is important to note that these
charges are not uniquely defined: if, for example, $\tau^{z}_{j}$ and
$\tau^{z}_{j+1}$ are conserved charges localized near lattice site $j$ with
$[H,\tau^{z}_{j}]=[H,\tau^{z}_{j+1}]=0$ then $\tau^{z}_{j}\tau^{z}_{j+1}$ is
also local and conserved, $[H,\tau^{z}_{j}\tau^{z}_{j+1}]=0$.

One way to obtain these charges is to try to perturbatively construct them in
the limit of strong disorder [40, 41, 43]. In this limit, one can think of
these charges as the local operators $\sigma^{z}_{i}$ smeared out by the
hopping and dressed by the interactions. As in our previous considerations,
there is then again an important difference between the non-interacting and
the interacting case: While in the former the one-particle operators
$\sigma^{z}_{i}=c^{\dagger}_{i}c_{i}-1$ remain one-particle operators, they
become many-particle operators in the latter case. As we will see, the
perturbative construction of local conserved charges will lead us again to
multiple commutators which are closely related to the ones studied in the
previous sections.

The Hamiltonian of a quantum system can be diagonalized by a unitary
transformation

$H=\exp(T)H\exp(-T)=\sum_{k=0}^{\infty}\frac{[T,H]^{(k)}}{k!}$ (IV.1)

with $T^{\dagger}=-T$ being an anti-unitary operator. For a Hamiltonian
$H=H_{0}+V$ where $H_{0}$ is diagonal in the operator basis chosen, we can
instead of a full diagonalization demand that the new Hamiltonian $H^{(1)}$ is
diagonal only to first order in the perturbation $V$. This is the Schrieffer-
Wolff transformation [46, 47]. Going one step further, we can consider
successive Schrieffer-Wolff transformations where in each step the
perturbation is eliminated to leading order. This is closely related to the
approach in Ref. [40, 41] which tries to generalize a proof of Anderson
localization to the many-body case. In these works, the author is mostly
concerned with the impact of resonances and splits the off-diagonal part of
the Hamiltonian into a resonant and a non-resonant part. We will discuss the
impact of resonances as well but will point out an even more fundamental issue
which calls into question the entire perturbative construction of local
conserved charges in the interacting case even for disorder configurations
where no resonances are present.

To be specific, we consider again the XXZ Hamiltonian (III.1). We note that in
Refs. [40, 41] the transverse field Ising model with random couplings and
random fields was considered instead. However, the arguments made in the
following are quite general and will apply to the latter model as well [48].
We want to investigate the case of large disorder where the magnetic field and
the interaction terms constitute $H_{0}$ with $[H_{0},\sigma^{z}_{i}]=0$ and
the hopping terms are the perturbation. The condition the operator $T$ has to
fulfill to lowest order is then given by

$V+[T,H_{0}]=0$ (IV.2)

and the new, transformed Hamiltonian is obtained as

$H^{(1)}=H_{0}+\frac{1}{2}[T,V]+\frac{1}{3}[T,[T,V]]+\cdots\,.$ (IV.3)

In the Schrieffer-Wolff transformation, one typically stops at the first order
correction to $H_{0}$. To justify this perturbative approach, we can rescale
the Hamiltonian (III.1) by $1/D$ which changes the energy scales but leaves
the physical properties unchanged. Then, the hopping terms have an amplitude
$\bar{J}^{(0)}_{j,j+1}=1/D$ and the interaction terms an amplitude
$\bar{\Delta}=\Delta/D$ which are both small parameters for large $D$. After
this rescaling, the random magnetic fields are always given by
$h_{j}\in[-1,1]$.

###  IV.1 Non-interacting case

First, we consider again the non-interacting case. An ansatz for the operator
$T$ can be obtained by commuting $H_{0}$ and $V$. This leads to
$T=\sum_{j}A_{j}(S^{+}_{j}S^{-}_{j+1}-S^{+}_{j+1}S^{-}_{j})$. The amplitudes
$A_{j}$ can then be determined from Eq. (IV.2) leading to

$T=\frac{1}{2}\sum_{j}\frac{\bar{J}^{(0)}_{j,j+1}}{h_{j}-h_{j+1}}(\sigma^{+}_{j}\sigma^{-}_{j+1}-\sigma^{+}_{j+1}\sigma^{-}_{j})\,.$
(IV.4)

We note that the Schrieffer-Wolff operator $T$ is thus proportional to the
spin current operator $\mathcal{J}=\sum_{i}j_{i}$ which can be obtained from
the lattice continuity equation [49]

$\partial_{t}\sigma^{z}_{i}=-i[\sigma^{z}_{i},H]=-(j_{i}-j_{i-1})\,.$ (IV.5)

We note, furthermore, that $\overline{|h_{j+1}-h_{j}|}=2/3$, i.e., the factor
$\bar{J}^{(0)}_{j,j+1}/|h_{j}-h_{j+1}|$ is small on average if the disorder
strength $D$ is large. However, if two neighboring sites have roughly the same
potential, then this factor can become large or might even diverge. Rare
resonances thus might destroy the perturbative construction if they
proliferate. We will come back to this point below.

After the Schrieffer-Wolff transformation, the new Hamiltonian is given by

$\displaystyle H^{(1)}$ $\displaystyle=$ $\displaystyle
H^{(0)}_{0}+\frac{1}{2}[T^{(0)},V^{(0)}]+\frac{1}{3}[T^{(0)},[T^{(0)},V^{(0)}]]$
(IV.6) $\displaystyle=$
$\displaystyle\sum_{j}\left\\{\bar{J}^{(1)}_{j,j+2}(\sigma^{+}_{j}\sigma^{z}_{j+1}\sigma^{-}_{j+2}+h.c.)+2h_{j}^{(1)}\sigma^{z}_{j}\right\\}$
$\displaystyle+$
$\displaystyle\sum_{j}\bar{J}^{(1)}_{j,j+1}(\sigma^{+}_{j}\sigma^{-}_{j+1}+h.c.)$

with $H^{(0)}_{0}=H_{0}$, $V^{(0)}=V$, and $T^{(0)}=T$. The coupling constants
and renormalized fields are given by

$\displaystyle\bar{J}^{(1)}_{j,j+2}$ $\displaystyle=$
$\displaystyle\frac{(\bar{J}^{(0)}_{j,j+1})^{2}}{2(h_{j+1}-h_{j})}+\frac{(\bar{J}^{(0)}_{j,j+1})^{2}}{2(h_{j+1}-h_{j+2})}$
(IV.7) $\displaystyle h_{j}^{(1)}$ $\displaystyle=$ $\displaystyle
h_{j}^{(0)}+\frac{(\bar{J}^{(0)}_{j,j+1})^{2}}{4(h_{j}-h_{j-k})}+\frac{(\bar{J}^{(0)}_{j,j+1})^{2}}{4(h_{j}-h_{j+1})}$

and $J^{(1)}_{j,j+1}\sim 1/D^{3}$ (the exact expression is lengthy and not
needed in the following). The first two terms correspond to those generated by
$[T^{(0)},V^{(0)}]$ while the renormalized nearest-neighbor hopping
$J^{(1)}_{j,j+1}$ stems from the next higher commutator
$[T^{(0)},[T^{(0)},V^{(0)}]]$. I.e., the transformation to lowest order
eliminates the nearest-neighbor hopping processes but generates a next-nearest
neighbor hopping $\bar{J}^{(1)}_{j,j+2}\sim 1/D^{2}$ (note that
$\sigma^{+}_{j}\sigma^{z}_{j+1}\sigma^{-}_{j+2}+h.c.=c^{\dagger}_{j}c_{j+2}+h.c.$
in the fermionic language) and renormalizes the local fields $h^{(0)}_{j}\to
h^{(1)}_{j}$.

If we stop at the leading order of the transformation, then we obtain a
Hamiltonian $H^{(1)}$ which looks like the original Hamiltonian $H^{(0)}$ but
with the nearest-neighbor hopping replaced by a weaker next-nearest neighbor
hopping. We can then again define a Schrieffer-Wolff transformation to
eliminate these next-nearest hopping terms which has the form

$T^{(k)}=\frac{1}{2}\sum_{j}\frac{\bar{J}^{(k)}_{j,j+k+1}}{h^{(k)}_{j+k+1}-h^{(k)}_{j}}(\sigma^{+}_{j}\sigma^{z}_{j+1}\dots\sigma^{z}_{j+k}\sigma^{-}_{j+k+1}-h.c.)$
(IV.8)

with $k=1$. This will again, to leading order, renormalize the local fields
and generate a third-nearest neighbor hopping which again can be eliminated to
leading order by a the Schrieffer-Wolff transformation (IV.8) with $k=2$.
I.e., we can define a series of successive Schrieffer-Wolff transformations
which lead to longer and longer range interactions which are, however,
suppressed as $D^{-2k}$. At the same time, they define an exponentially
converging series for the renormalized random fields $h^{(k)}_{j}$. If there
are no resonances, then we can drop the hopping term when all the amplitudes
$J^{(k)}$ have fallen below some threshold, giving us an effective Hamiltonian
with renormalized magnetic fields. In this approximation, which can be
justified at very strong disorder and which corresponds to summing up exactly
certain classes of diagrams, the conserved local charges are just the
$\sigma^{z}_{j}$ themselves but with renormalized amplitudes. This can be seen
as a consistency check that at very large disorder such an effective
description exists because the hopping is an irrelevant perturbation.

In the non-interacting case, it is also well understood that distance
resonances have a vanishingly small probability to proliferate [1, 3, 4]. In
our construction, this scaling argument arises as follows: Two lattice sites
$j$ and $j+\ell$ a distance $\ell$ apart are coupled by an effective long-
range hopping $\bar{J}^{(\ell-1)}_{j,j+\ell}(c_{j}^{\dagger}c_{j+\ell}+h.c.)$
with an effective amplitude obtained by $\ell-1$ successive Schrieffer-Wolff
transformations. The amplitude of this long-range effective hopping term is
given by $\bar{J}^{\ell-1}\sim D^{-2\ell+2}\times
f[\\{(h^{\ell-1}_{j}-h^{\ell-1}_{k})^{-1}\\}]$ where $f$ is some function of
the inverse differences of the random fields within the cluster of $\ell$
sites. Now the smallest pairwise difference on this cluster on average scales
as $\ell^{-2}$. Thus the effective long-range hopping amplitude scales as
$\bar{J}^{(\ell-1)}_{j,j+\ell}\sim D^{-2\ell}\ell^{2}$ and is thus typically
exponentially suppressed. Therefore resonances have a vanishingly small
probability to proliferate and to destroy the perturbative construction above.
This is also supported by our symbolic calculations on finite chains of length
$N$. We do successive lowest order Schrieffer-Wolff transformations and stop
the transformations if the total 1-norm of the off-diagonal part is below some
threshold value $\varepsilon$. We have $\sim N$ possible off-diagonal terms of
a given range and if the amplitudes of all those terms, on average, are
suppressed as $D^{-2k}$ at order $k$ of the transformation then our condition
to stop is $ND^{-2k}<\varepsilon$. We thus expect that we need
$k\sim\ln(N/\varepsilon)$ many transformation steps. This is consistent with
our observations.

Quite generally we may ask why we expect a perturbative construction of the
local conserved charges to work in the non-interacting case? In the lowest-
order construction discussed above, the Hamiltonian $H^{(k)}$ always has $3N$
terms for a chain of length $N$ for every $k$. I.e., the number of terms does
not increase during successive Schrieffer-Wolff transformations. If we take
higher order commutators in the transformation into account, then the number
of terms will increase, however, the Hamiltonian $H^{(k)}$ and the
transformation $T^{(k)}$ will always consist of a sum of one-body operators
and, for a chain of length $N$, there are only $N^{2}$ distinct one-body
operators. Therefore the problem of finding local conserved charges is at most
quadratic, and, in any finite order approximation, one in fact always deals
with a constant, non-increasing number of terms in the transformed
Hamiltonians.

###  IV.2 Interacting case

The interacting case, on the other hand, is completely different and none of
the above arguments hold. We are, however, able to generalize Eq. (IV.4) and
to obtain the exact Schrieffer-Wolff transformation for the XXZ chain and
arbitrary interaction strength $\Delta$. We find

$\displaystyle T$ $\displaystyle=$
$\displaystyle\sum_{j}\frac{\bar{J}^{(0)}_{j,j+1}(\sigma^{+}_{j}\sigma^{-}_{j+1}-h.c.)}{2(h_{j}-h_{j+1})(2\bar{\Delta}^{2}-2(h_{j}-h_{j+1})^{2})}$
(IV.9) $\displaystyle\times$
$\displaystyle[\bar{\Delta}^{2}-2(h_{j}-h_{j+1})^{2}+\bar{\Delta}^{2}\sigma^{z}_{j-1}\sigma^{z}_{j+2}$
$\displaystyle+\bar{\Delta}(h_{j}-h_{j+1})(\sigma^{z}_{j-1}-\sigma^{z}_{j+2})]$

with $\bar{J}^{(0)}_{j,j+1}=1/D$ and $\bar{\Delta}=\Delta/D$. We note that for
$\Delta\neq 0$, the transformation $T$ consists of $8N$ distinct operators for
a chain of length $N$, including operators acting on $2,3$ and $4$ neighboring
sites.

If we stop the transformation (IV.3) at first order, the transformed
Hamiltonian $H^{(1)}$ already has $29N$ distinct operators if the chain length
$N$ is sufficiently large. This includes terms acting on $5$ neighboring sites
but also terms such as $\sigma^{z}_{j}\sigma^{z}_{j+1}$ and
$\sigma^{z}_{j}\sigma^{z}_{j+2}\sigma^{z}_{j+3}$. Even for very strong
disorder, the transformation thus does not just lead to a renormalization of
the random fields but necessarily leads to a broadening of the diagonal basis.
Symbolically, we can also construct $T^{(1)}$ to eliminate the off-diagonal
terms in $H^{(1)}$. We find that $T^{(1)}$ already consists of $786N$ distinct
operators. The transformed Hamiltonian $H^{(2)}$ then has $30618N$ many terms.
I.e., even if we stop the Schrieffer-Wolff transformation at lowest order, the
number of terms in $T^{(k)}$ and $H^{(k)}$ is growing exponentially.
Furthermore, there is no distinct set of diagrams which we can simply sum up
to all orders. Exponentially many new diagrams are created at each order.

The related and most fundamental problem, however, is that lattice sites $j$
and $j+\ell$ a distance $\ell$ apart will be coupled by $\sim 3^{2}4^{\ell-1}$
distinct operators because on each site we can put in principle one of the
operators $\\{\sigma^{x},\sigma^{y},\sigma^{z},\mathbb{I}\\}$ but the string
of operators of length $\ell+1$ has to start and end with one of the three
Pauli matrices. While not all of these combinations of operators will be
realized in a microscopic model such as the XXX chain, there will be
exponentially many distinct operators which couple these two sites, see Fig.
15. Thus there is, in general, no reason to expect this perturbative
construction to work because the overall effective amplitude of terms coupling
these two sites will be $\bar{J}_{j,j+\ell}\sim
D^{-\ell}4^{\ell}\sim\mathcal{O}(1)$. I.e., there is no small parameter
anymore. This reasoning is supported by the results of symbolic calculations
shown in Fig. 15. Here we have calculated the sum $s_{\ell}^{(1)}(k)$ of the
amplitudes of all terms obtained in the Schrieffer-Wolff transformation (IV.9)
to order $k$ which have support on $\ell$ sites, i.e., which connect sites a
distance $\ell-1$ apart. For a periodic chain with even length $N$, there are
$\sim N$ distinct pairs of sites which are a distance $1,\cdots,N/2$ apart.
I.e., for a specific pair the overall coupling is on average $\sim
s_{\ell}^{(1)}(k)/N$ and, as can be seen in Fig. 15, is in general not
exponentially decreasing with $\ell$. We note that the system considered in
Fig. 15 has only $N=20$ sites so that only results up to $\ell\sim 10$ are
relevant. That the number of terms generated in such a transformation with
support on $\ell$ sites is exponentially increasing with $\ell$ makes it very
questionable that a construction of local conserved charges in this
perturbative manner is a controlled approximation. This is a major issue even
before thinking about the problem of resonances and thermal inclusions
(regions for a given disorder configuration which, by chance, are only weakly
disordered). It seems to us that the exponential increase of the number of
terms in such transformations is a point which has not been properly addressed
in the proof of MBL in Refs. [40, 41] where the author seems to be only
concerned about the effect of resonances. While it is not impossible that for
disorder configurations without resonances and thermal inclusions the
perturbative construction of local charges is working in the thermodynamic
limit—for example, because of interference effects between different terms
coupling the same sites—any proof needs to provide rigorous arguments that the
exponentially many terms do not overwhelm the exponential small parameter in
the expansion. Such arguments seem to be missing so far. We note that for a
finite system the construction will always eventually converge because the
finite length provides a cutoff and the norm $s_{\ell}^{(1)}(k)$ will
eventually drop off exponentially with $\ell$ as can also be seen in Fig. 15.

The destabilizing effects of resonances and thermal inclusions have received
considerable attention in the literature already [50, 51, 35]. From the
perspective of perturbative Schrieffer-Wolff transformations, they lead to the
following additional problems: The renormalized couplings for an operator
acting on $\ell$ lattice sites will in its denominator typically contain a
linear combination of all the random fields on these $\ell$ sites. Given that
there are exponentially many such operators, the simple scaling argument for
the non-proliferation of resonances in the Anderson case is clearly no longer
applicable in the interacting case. In particular, thermal inclusions—regions
of the lattice which have similar potentials—lead to a further destabilization
of the perturbative transformation. This can be seen, for example, in symbolic
calculations of the unitary transformation (IV.1) with the operator (IV.9) to
higher orders. In configurations where such regions exist, the transformation
generates long-range terms with large amplitudes even if the disorder $D$ is
extremely strong. If, on the other hand, we exclude configurations with such
regions and with any pairwise resonances by hand, then the amplitudes of any
individual operator typically decay with the length of its support. Another
manifestation of this issue can be seen when performing $k$ successive lowest
order Schrieffer-Wolff transformations for a finite chain of length $N$. We
want to stop the transformations again once the 1-norm of the off-diagonal
part is below some threshold $\varepsilon$. In contrast to the non-interacting
case, we now have exponentially many terms connecting sites a given distance
apart. If we assume that individual amplitudes in the off-diagonal part are
suppressed as $D^{-2k}$, then the condition for stopping the transformations
becomes $\text{e}^{N}D^{-2k}<\varepsilon$ and we would thus expect that we
need $k\sim N$ steps for $N$ sufficiently large. This is, however, not what we
observe for typical disorder configurations. The results of our symbolic
calculations point instead to an increase of the number of transformations
$k\sim\text{e}^{N}$. This implies that some amplitudes in the off-diagonal
part are $\sim\mathcal{O}(1)$ and that we have to eliminate almost all of the
off-diagonal terms to achieve convergence.

In summary, the iterative construction of local conserved charges cannot
simply be transferred from the non-interacting to the interacting case. There
are two main issues: Even if one considers only a lowest order transformation,
one is not dealing with a constant number of operators but rather with an
exponentially increasing number of operators, both in the Schrieffer-Wolff
transformation $T^{(k)}$ as well as in the transformed Hamiltonian $H^{(k)}$.
Even if one tries to only perform a single unitary transformation (IV.1) with
the operator $T$, the number of distinct operators in the transformed
Hamiltonian is increasing exponentially with the support of the operator which
thus can overwhelm the exponential decay of the amplitude of individual
operators with the length of their support. We note that different but related
arguments for the failure of a perturbative construction of conserved charges
in the interacting case have recently been made in Ref. [52]. Furthermore,
regions with similar lattice potentials result in individual operators with
large support which, nevertheless, have large amplitudes even for very strong
disorder. The latter issue is related to the well-known avalanche instability
[50, 51, 53, 35].

##  V Conclusions

We have considered operator growth in one-dimensional lattice models and its
connection to localization from two different perspectives. First, we studied
the multiple commutator $[H,A]^{(k)}$, which appears in the Euclidean time
evolution, for a local operator $A$ and a Hamiltonian $H$ which is a sum of
nearest-neighbor terms. For a generic system, we have shown that the number of
distinct terms generated is exponentially increasing with the order $k$ of the
commutator. This leads to a strict bound for the commutator norm $s(k)$ which
grows faster than exponential with $k$, a result which was first demonstrated
in Ref. [30]. Since this result is based on a simple counting of connected
clusters generated by the commutator, it is expected that this bound is
asymptotically tight for a generic system. We then provided strict bounds for
two special cases. First, we considered the Hamiltonian $H$ of a free-fermion
system. In this case, all the terms generated by the commutator of the
Hamiltonian with a one-body operator $A$ will also always be one-body
operators. For a short-range Hamiltonian, the spatial range of operators
generated at order $k$ will be proportional to $k$. The number of distinct
terms in the commutator can therefore only grow $\sim k^{2}$ as compared to
the exponential growth in the general, interacting case. We also proved an
exponential bound for the commutator norm in this case. Second, we considered
a localized system, either interacting or non-interacting, and showed that an
exponential bound for $s(k)$ exists in this case as well. Physically, this can
be understood as follows: a localized system has a set of local conserved
charges and the Hamiltonian can be written in terms of these charges.
Commuting the Hamiltonian with a local operator $A$ then only involves those
charges localized near the site the operator $A$ is acting on. Thus, up to
exponential tails, only operators with support on a number of sites which is
less than the localization length will significantly contribute to the norm.

To distinguish a localized interacting or non-interacting system from a non-
localized free fermion system, it is important to look beyond the overall norm
which has an exponential bound in both cases. Here we have shown that for a
non-localized free-fermion system at large commutator order $k$, the norm
$s_{l}(k)$ of terms which have support on $l$ sites will grow exponentially
with $k$ and is independent of the length of the support $l$. In contrast, the
norm scales as $s_{l}(k)\sim\exp(k)\exp(-l)$ for a localized system. I.e., at
a fixed order $k$, the contribution of terms in the commutator is
exponentially suppressed with the length of their support $l$.

We thus have proven sharp criteria to distinguish between free-fermionic
localized and non-localized and interacting localized and non-localized
systems based on the growth of the total commutator norm $s(k)$ and the
contributions $s_{l}(k)$ to this norm by operators with support on $l$ sites.
To illustrate our results, we then considered specific microscopic examples,
namely the XX chain (nearest-neighbor free-fermion model), the Anderson model,
the Aubry-André model, and the XXX (Heisenberg) chain with and without
quenched disorder. For the XX case we have demonstrated that a complete, exact
result for the commutator $[H,\sigma^{z}_{0}]^{(k)}$ to all orders $k$ can be
derived. These exact results confirmed and illustrated our general results
derived earlier. For the Anderson model, we derived an exponential bound for
the total norm which does depend on the disorder strength $D$. Based on a
graph theoretical approach, we also showed that $s_{l}(k)\sim\exp(-l)$ for
fixed $k$. I.e., the contribution of terms in the commutator to the total norm
is indeed exponentially suppressed with the length $l$ of their support. As
another example, we considered the Aubry-André model which has a transition
from a non-localized to a localized phase at some finite strength of the
applied incommensurate lattice potential. Importantly, we could demonstrate
that the non-localized and localized phases of the model can clearly be
distinguished by considering the norm $s_{l}(k)$.

Having checked our results for non-interacting cases, we then turned to the
Heisenberg (XXX) chain with quenched disorder. This model has been suggested
to have a transition from an ergodic to a many-body localized phase at some
finite disorder strength $D$. In particular, for the MBL phase an effective
Hamiltonian written in terms of localized conserved charges is supposed to
exist as in the Anderson case but with couplings between these charges which
are exponentially decaying with distance. If true, then our proof for
localized phases applies. Thus $s(k)$ can grow at most exponentially in such a
phase and $s_{l}(k)\sim\exp(-l)$ at fixed $k$. We used symbolic calculations
of the commutator and were able to reach order $k=16$ at which point we are
considering terms located on a lattice of $2k+1=33$ sites. This is
substantially larger than what can be done using exact diagonalizations. We
found that the growth of the norm of the commutator $s(k)$ appears to be
faster than exponential and consistent with the generic bound for non-
localized systems for all considered disorder strengths up to $D=25$.
Similarly, we find no indication that $s_{l}(k)\sim\exp(-l)$. To the contrary,
we find that $s_{l}(k)/s_{1}(k)$ has a maximum which appears to shift to
larger and larger $l$ with increasing $k$. At the very least, we claim that
our results clearly demonstrate that there is no convincing case for an MBL
phase for system sizes up to $33$ sites and disorder strengths up to $D=25$.

As a second approach, we considered unitary Schrieffer-Wolff transformations
of the entire Hamiltonian to perturbatively eliminate the hopping terms and
construct local conserved charges. We showed that such an approach will
succeed in the Anderson case because the Hamiltonian in this case consists of
one-body operators and this property is not changed by unitary
transformations. Consequently, there are only two terms possible in the
effective Hamiltonian which can couple two sites $j$ and $k$:
$c^{\dagger}_{j}c_{k}$ and $c_{j}c_{k}^{\dagger}$. This allows to define
consecutive lowest order Schrieffer-Wolff transformations and to sum up the
entire series. In this approximation, valid at very strong disorder, the local
conserved charges are simply the $\sigma^{z}_{j}$ operators with renormalized
fields. The longer-range hoppings generated by the transformations have
amplitudes which decay exponentially with the range of the hopping process.
The minimum gap between random fields on a cluster of length $\ell$ on the
other hand—which enters the effective hopping in the denominator—scales as
$\ell^{-2}$ and thus cannot overcome the general suppression of the hopping
process which scales as $D^{-2\ell}$. This means, in particular, that
resonances cannot proliferate.

All of this changes drastically once interactions are included. In particular,
there are then $\sim 4^{\ell}$ different terms which can connect sites a
distance $\ell$ apart. We confirm this exponential growth with the length of
the support by symbolic calculations. A priori, there is thus no reason to
believe that a perturbative construction of the local conserved charges by
unitary transformations is still possible in the interacting case. Even if we
assume that there are no resonances or thermal inclusions (regions with random
fields of similar strengths), there is a fundamental problem with this
construction: while the amplitude of an individual term connecting sites a
distance $\ell$ apart will still be $\sim D^{-2\ell}$, there are exponentially
many such terms so that the total effective amplitude is $\sim
4^{\ell}D^{-2\ell}\sim\mathcal{O}(1)$. The small parameter in the perturbative
construction is lost. Our symbolic calculations do confirm this issue and
indeed show that the total effective amplitude is not decaying with the
support of the operator in contrast to the non-interacting case. As far as we
can tell, this issue is not addressed at all in the proof of MBL by Imbrie
[40, 41] which uses the same kind of unitary transformations but seems to be
only concerned about resonances. The latter is also a much more severe problem
in the interacting case. Even if we assume that the perturbative construction
is convergent for carefully chosen disorder configurations without resonances
and thermal inclusions due to some miraculous cancellation of contributions
stemming from different commutator paths, many-body resonances will destroy
the approach when present. A long-range term typically involves operators on
all sites between the left most and the right most site. Consequently, the
renormalization will involve linear combinations of all the random fields on
this cluster and there are exponentially many of such terms. Clearly, the
argument for the non-proliferation of resonances made in the non-interacting
case is then no longer valid. At the very least, we believe that we have shown
that there is no reason to believe that the proof in Refs. [40, 41] is valid.
It seems to ignore that, in general, the number of terms connecting distant
sites grows exponential with the distance between the sites. No argument is
made why the construction remains controlled by a small parameter despite the
exponential number of terms and why resonances do not proliferate. Our
symbolic calculations for finite chains and finite orders of the
transformation suggest that such an argument, in fact, cannot be made and that
a perturbative construction of local conserved charges in this way is
impossible except for $D=\infty$ where the model is trivially localized.

While we have concentrated here on the Heisenberg chain with quenched
disorder, we note that in the attempted proof of MBL by Imbrie [40, 41] the
transverse Ising chain with random couplings and fields in considered instead.
For the latter model, it has recently been rigorously proven [48] that the
norm of a local operator when commuted $k$-times with the Hamiltonian grows
faster than exponential with $k$ and, in fact, asymptotically like the general
bound (II.12). Here, we have shown that if a local operator remains quasi-
local when commuted with a Hamiltonian then the total norm cannot grow faster
than exponential. If we negate this statement, then this means that if the
commutator norm of a local operator with a Hamiltonian grows faster than
exponential — as in the transverse Ising case — then the commutator is
necessarily not a quasi-local operator. It can at best be algebraically
localized. I.e., our results come close to proving that an MBL phase does not
exist in the transverse Ising chain. The only assumption we have made is that
for the commutator to remain quasi-local, contributions to the connected
clusters coming from sites far away from the local operator are exponentially
suppressed, see App. A. This seems to be a physical requirement for a
localized model. Mathematically it is also hard to see how a local operator
can otherwise remain quasi-local. This would then require some miraculous
cancellation between different paths which has to happen for any local
operator and most disorder configurations. Making these arguments even more
rigorous might be the last remaining step in proving that localization is not
possible in the many-body case.

###### Acknowledgements.

J.S. acknowledges discussions with D. Sels and M. Kiefer-Emmanouilidis and
funding by the NSERC Discovery program and by the DFG via the research unit
FOR 2316.

##  Appendix A Localized Models

We can strengthen our argument from section II.3 so that rather than a hard
cutoff for contributing connected clusters, we consider contributions from
bonds increasingly distant from our original site to decay exponentially. We
wish to insert this condition into the counting of connected clusters. To do
this, first note that the Stirling number of the second kind can be written as

$S(k,j)=\sum_{k_{0}+\ldots+k_{j-1}=k}\frac{1}{j!}\prod_{i}\frac{1}{a_{i}!}{{k}\choose{k_{0},k_{1},\ldots,k_{j-1}}}$
(A.1)

where each $k_{i}$ is an integer greater than or equal to 1. The term in
brackets is the multinomial coefficient which is defined as

${{k}\choose{k_{0},k_{1},\ldots,k_{j-1}}}=\frac{k!}{k_{0}!k_{1}!\cdots
k_{j-1}!}.$ (A.2)

Eq. (A.2) can be interpreted as the number of ways of distributing $k$ objects
into $j$ boxes with $k_{i}$ objects in box $i$, $0\leq i\leq j-1$. Unlike the
Stirling number of the second kind $S(k,j)$, we are now specifying the number
of objects in each of the $j$ boxes, which is why we have the sum in Eq.
(A.1). Further, we now account for the order of the boxes; eliminating this
constraint leads to the factor of $\frac{1}{j!}$ in Eq. (A.1). Finally, the
product of $\frac{1}{a_{i}!}$ factors corrects for the fact that the
multinomial coefficient treats repeated integers in the set
$\\{k_{1},\ldots,k_{j}\\}$ as separate cases while $S(k,j)$ does not. To
better illustrate this formula, we provide an example of its usage.

Example: Consider $k=6$, $j=4$. Then our solutions to $k_{1}+\ldots+k_{j}=k$
will involve permutations of $\\{k_{1},k_{2},k_{3},k_{4}\\}=\\{1,1,1,3\\}$ or
$\\{1,1,2,2\\}$. There will be $4!=24$ terms in the sum for each of these
possibilities, but since the multinomial coefficient does not change value
when swapping any of the $k_{i}$’s, the $\frac{1}{j!}$ factor turns this into
a sum of the two multinomial coefficients with
$\\{k_{1},k_{2},k_{3},k_{4}\\}=\\{1,1,1,3\\},\\{1,1,2,2\\}$. However, we still
have the correcting $\frac{1}{a_{i}!}$ product. There are three repeated
integers in the set $\\{1,1,1,3\\}$ so we take
$\prod_{i}\frac{1}{a_{i}!}=\frac{1}{3!}$, while there are two sets of two
repeated integers in the set $\\{1,1,2,2\\}$ so we take
$\prod_{i}\frac{1}{a_{i}!}=\frac{1}{2!2!}$. The r.h.s. of Eq. (A.1) then
becomes

$\frac{1}{3!}{{6}\choose{1,1,1,3}}+\frac{1}{2!2!}{{6}\choose{1,1,2,2}}=65$
(A.3)

and in fact, $S(6,4)=65$.

We now denote by $\tilde{S}(k,j)$ our ‘adjusted’ Stirling number which inserts
the requirement that contributions from sites of increasing distance from our
initial site are exponentially decaying. Let $k_{n}$ denote the number of
bonds at a distance $n$ from our initial site, where a distance of zero
corresponds to a bond involving our initial site. We want exponential decay of
the ratio of the contributions from bonds at $k_{n}$ to bonds at $k_{0}$. To
insert this condition, we multiply by a factor of $e^{-\xi nk_{n}}$ for each
$k_{n}$ and for any constant $\xi>0$. We note that this is normalized so that
at $k_{0}$ this factor is 1, meaning the initial site is not reduced. This
separate treatment of the $k_{n}$’s, which represent the number of times each
bond appears in the connected cluster, is why the form (A.1) of $S(k,j)$ is
desirable. Applying the decaying factors, we get

$\displaystyle\tilde{S}(k,j)$ (A.4) $\displaystyle=$
$\displaystyle\sum_{k_{0}+\ldots+k_{j-1}=k}\frac{1}{j!}\prod_{i}\frac{1}{a_{i}!}{{k}\choose{k_{0},\ldots,k_{j-1}}}\prod_{n=0}^{j-1}e^{-\xi
nk_{n}}$ $\displaystyle\leq$
$\displaystyle\sum_{k_{0}+\ldots+k_{j-1}=k}{{k}\choose{k_{0},\ldots,k_{j-1}}}\prod_{n=0}^{j-1}e^{-\xi
nk_{n}}$ $\displaystyle=$ $\displaystyle\left(\frac{e^{\xi}(1-e^{-\xi
j})}{e^{\xi}-1}\right)^{k}.$

We also loosened the condition on the $k_{i}$’s in the second line so that
they are now greater than or equal to $0$ instead of 1 so that we may apply
the multinomial theorem in the third line. Since this may only increase the
sum due to allowing more solutions to $k_{0}+\ldots+k_{j-1}=k$, the upper
bound is still valid. Then the adjusted number of connected clusters, with the
inserted requirement of the exponentially decreasing contribution with cluster
length, is bounded from above by

$\displaystyle\sum_{\\{I_{1},\cdots,I_{k}\\}}\\!\\!\\!\\!\\!\\!\\!{\vphantom{\sum}}^{\prime}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{k}2^{j}\tilde{S}(k,j)$ (A.5)
$\displaystyle\leq$ $\displaystyle 2(1-e^{-\xi
k})^{k}\left(\frac{e^{\xi}}{e^{\xi}-1}\right)^{k}(2^{k}-1)$
$\displaystyle\approx$ $\displaystyle
2\left(\frac{2e^{\xi}}{e^{\xi}-1}\right)^{k}$

where in the last line we took the large $k$ limit. Defining the constant
$\alpha=\frac{2e^{\xi}}{e^{\xi}-1}$, we therefore obtain the bound

$||[H,A]^{(k)}||\leq 2||A||(2\alpha\tilde{J})^{k}\,.$ (A.6)

We note that this more careful calculation is qualitatively fully consistent
with Eq. (II.21) in the main text where we made the simplifying assumption
that clusters longer than the localization length do not contribute at all to
the norm. We thus conclude that the norm growth of the commutator in any
localized model is at most exponential.

##  Appendix B Graphs

The growth of the 1-norm of the commutator
$\left[H,\sigma_{0}^{z}\right]^{(k)}$ can also be understood using graphs both
without and with disorder. Here we provide some more details how such graphs
can be constructed and how they have to be read.

###  B.1 No disorder

We use a rooted, directed, weighted, infinite graph to model the commutator
growth. In the graph, each node represents classes of terms formed by the
commutator and each directed edge represents connections between classes of
terms formed via commutation with the Hamiltonian. For example, if commuting a
term of class $A$ with $H$ produces a term of class $B$ then an arrow is drawn
on the graph from $A$ to $B$. Further, the 1-norm of terms of class $B$ in
that commutation is assigned to the weight of that edge. The classes are
organized such that each term of that class produces the same result when
commuting with $H$.

In the non-interacting, non-disordered case, whose exact commutator structure
is given by Eq. (III.1.1) and Eq. (III.1.1), we find that the growth is
modelled by the left column of the graph of Fig. 4. By convention, we label
the classes by the string of Pauli operators beginning with $\sigma^{x}$, and
we omit index labels with the understanding that a string of operators acts on
consecutive sites. Note that terms are invariant with respect to commutation
with $H$ under the transformation $\sigma^{x}\leftrightarrow\sigma^{y}$, so we
group terms of the form $\sigma^{x}\cdots\sigma^{y}$ and
$\sigma^{y}\cdots\sigma^{x}$ into the same class, and
$\sigma^{x}\cdots\sigma^{x}$ and $\sigma^{y}\cdots\sigma^{y}$, but not
$\sigma^{x}\cdots\sigma^{y}$ and $\sigma^{x}\cdots\sigma^{x}$.

Enumerative calculations on the graph Fig. 4 now allow us to calculate the
1-norm growth and spatial structure of the commutator
$\left[H,\sigma_{0}^{z}\right]^{(k)}$. We always begin with a single
$\sigma^{z}$ operator which we label node 1. We then commute with the
Hamiltonian $k$ times which corresponds to a walk of length $k$ on the graph
starting at node 1. To get the total 1-norm, we must sum over all weighted
walks of length $k$ which begin at node 1 with no restrictions on the end
node. By a weighted walk we mean a walk through the edges of the graph in
which the running count is multiplied by the weight of the traversed edge.

It is simple to count the total number of weighted walks on the graph Fig. 4
for $D=0$ (left column) if we make the following observation. If on the $i$-th
step we are at node 1, then we have a single option for the next node, and
this traversal has weight 8. Thus, the weighted number of options is $8\cdot
1=8$. If we are at any other node, there are two options for the next node,
and each traversal has weight 4, so the weighted number of options is $4\cdot
2=8$. Since the number of weighted options is thus equal at each step and
there are $k$ total steps, the 1-norm must be equal to $8^{k}$, reproducing
our result of Eq. (III.8).

We label node $l$ to be the node corresponding to the class of terms which are
a string of $l$ Pauli matrices. The support on $l$ sites then corresponds to
the weighted number of walks which begin at node 1 and end at node $l$. We
first consider the support on one site. We note that any walk that begins and
ends at node 1 on the graph Fig. 4 must have even $k$ and traces out a Dyck
path. Thus, the total unweighted number of walks is equal to
$C\left(\frac{k}{2}\right)$, where $C(n)$ in this case is the Catalan number.
However, due to the edge leaving the first node having a different weight, to
find the total weighted number of Dyck paths one must use the following
formula

$s^{(1)}_{1}(k)=\sum_{m}w(m,k)p(m,k)$ (B.1)

where $s^{(1)}_{1}(k)$ means the 1-norm of terms with support on one site as a
function of $k$, $w(m,k)$ is the weight of a path as a function of $m$ and $k$
where $m$ is the number of traversals of the single edge of weight 8, and
$p(m,k)$ is the number of paths as a function of $m$ and $k$. Note that we use
the term paths instead of walks here to refer to the fact that these are Dyck
paths. It is clear that the bounds on $m$ are $1\leq m\leq\frac{k}{2}$ and
that $w(m,k)=4^{k-m}8^{m}=4^{k}2^{m}$. We calculate $p(m,k)$ by partitioning
the Dyck path into $m$ parts according to intersections with the $x$-axis. We
then have

$\displaystyle p(m,k)$ $\displaystyle=$
$\displaystyle\sum_{i_{1}+\ldots+i_{m}=\frac{k}{2}}\;\prod_{j=1}^{m}C(i_{j}-1)$
(B.2) $\displaystyle=$ $\displaystyle\frac{m}{k-m}{{k-m}\choose{\frac{k}{2}}}$

where we used Catalan’s k-fold convolution formula. This leads to

$s^{(1)}_{1}(k)=4^{k}\sum_{m=1}^{\frac{k}{2}}2^{m}\frac{m}{k-m}{{k-m}\choose{\frac{k}{2}}}=4^{k}{{k}\choose{\frac{k}{2}}}.$
(B.3)

The support on 2 sites is then easily found using a recursive strategy as

$s^{(1)}_{2}(k)=\frac{1}{4}s_{1}(k+1)=2^{2k+1}{{k}\choose{\frac{k-1}{2}}}\,.$
(B.4)

We see that this reproduces our previous result Eq. (III.1.1) with $l=2$.
Another recursive equation gives the support on 3 sites

$\displaystyle s^{(1)}_{3}(k)$ $\displaystyle=$
$\displaystyle\frac{1}{4}s_{2}(k+1)-2s_{1}(k)$ $\displaystyle=$ $\displaystyle
2^{2k+1}\left[{{k+1}\choose{\frac{k}{2}}}-{{k}\choose{\frac{k}{2}}}\right]=2^{2k+1}{{k}\choose{\frac{k-2}{2}}}.$

We see that this also agrees with our previous result Eq. (III.1.1) for $l=3$.
Finally, we get the following recursive equation for $l\geq 4$:

$s^{(1)}_{l}(k)=\frac{1}{4}s_{l-1}(k+1)-s_{l-2}(k).$ (B.6)

Simple inductive arguments show that the solution to Eq. (B.6) with base cases
(B.4) and (B.1) is indeed our previous result (III.1.1).

###  B.2 Disordered case

In the disordered case, we find that we have new classes of terms in the
commutator with edges given by the magnetic field values $h_{j}$. For the
Anderson case, we know that $h_{j}\in[-D,D]$, so the maximum absolute value of
the disorder at each site is $D$. For the Aubry-André model, the same
additional terms in the commutator appear but in this case the maximum value
is $\lambda/2$. This allows us to create a modified graph which contains an
additional column of terms, see Fig. 4, and which allows us to derive again an
upper bound of the 1-norm. We discuss the Anderson case here and note that the
Aubry-André case follows immediately with the replacement $D\to\lambda/2$.

Since $D$ is an upper bound on the disorder at each site, we can only compute
an upper bound on the 1-norm growth of the commutator. Following the same
strategy as for the non-disordered case, we see that at a given step, if we
are at node 1, then we have a weighted number of options of 8. If we are at
the $\sigma^{x}\sigma^{x}$ node, we have a weighted number of options of
$4+8D$. If we are at any other node, we have a weighted number of options of
$8+8D$. The maximum of these is $8+8D$, so we get the following upper bound on
the 1-norm growth of the commutator, keeping in mind that the first step must
have a contribution of 8:

$||[H,\sigma^{z}_{0}]||_{1}\leq 8^{k}(D+1)^{k-1}.$ (B.7)

This is the result Eq. (III.14). The bound is exponential as expected, with an
increasing disorder strength corresponding to an increased growth rate. We
also see that setting $D=0$ gives us $8^{k}$, the exact result in the case
without disorder.

Now we consider the spatial structure of the commutator. We can get a
qualitative sense of why localization occurs in this case. Recall that the
support on $l$ sites in the graphical language is the weighted number of walks
of length $k$ that start at node 1 and ends at node $l$. First note that in
this case there are two nodes for each $l$. Additionally, with the new terms
due to disorder, there will be more possible walks of length $k$ as one may
also traverse the edges with weight $8D$. However, traversing these edges does
not result in a node with an increased $l$. Therefore, in contrast to the non-
disordered case, there will be a greater number of walks that end at nodes
with smaller $l$, and the relative number of walks should decrease as $l$
increases. Furthermore, for $D\gg 1$ traversing the edges with weight $8D$,
which does not increase $l$, leads to a much larger contribution to the norm
than traversing those edges which connect nodes with different $l$’s. This
means that the $1$-norm of term with support on $l$ sites decays with $l$,
indicating localization.

For a more quantitative analysis, let us consider the upper bound on the
support at $l$ sites $s^{(1)}_{l}(k)$. Numerical calculations using weighted
adjacency matrices led to the following result for the support on one site

$\displaystyle s^{(1)}_{1}(k)$ $\displaystyle\leq$ $\displaystyle
2^{3k-1}D^{k-2}\sum_{i=0}^{\frac{k-2}{2}}C(i){{k-1}\choose{2i}}(2D)^{-2i}$
$\displaystyle=$ $\displaystyle
2^{3k-1}D^{k-2}{}_{2}F_{1}\left(\frac{1-k}{2},\frac{2-k}{2};2;\frac{1}{D^{2}}\right)$

where ${}_{2}F_{1}$ is the hypergeometric function. Note that this result is
derived for even $k$. For $k$ odd the support will be zero; one can infer this
from Fig. 4 since there is no possible walk on the graph beginning and ending
at the top node with an odd number of steps. The upper bound in Eq. (B.2) only
arises from the fact that $D$ is an upper bound on the local magnetic fields;
the enumerative calculations are exact. Similar numerical calculations show
that the support on $l>1$ sites is bounded from above by a hypergeometric
function as well but with a different leading term

$\displaystyle s^{(1)}_{l}(k)$ $\displaystyle\leq$ $\displaystyle
2^{3k+2-l}{{k-1}\choose{l-2}}D^{k-l+1}$ $\displaystyle\times$
$\displaystyle{}_{2}F_{1}\left(\frac{1-k}{2},\frac{2-k}{2};2;\frac{1}{D^{2}}\right).$

The upper bound is looser in this case because the enumerative calculations
are no longer exact. We are then interested in the ratio
$s^{(1)}_{l}(k)/s^{(1)}_{1}(k)$ in the limit $k\gg l\gg 1$. Since the upper
bound given for $s^{(1)}_{l}(k)$ is in fact looser than that given for
$s^{(1)}_{1}(k)$, we may bound the ratio by the ratio of the two given upper
bounds, leading to the result

$\displaystyle\frac{s^{(1)}_{l}(k)}{s^{(1)}_{1}(k)}$ $\displaystyle\leq$
$\displaystyle{{k-1}\choose{l-2}}(2D)^{3-l}$ (B.10) $\displaystyle\approx$
$\displaystyle\frac{1}{\sqrt{2\pi
l}}\left(\frac{ke}{l}\right)^{l}e^{-\frac{l^{2}}{2k}}(2D)^{3-l}$

This coincides with Eq. (III.17). As mentioned in the main text, there is a
maximum in this ratio which shifts to larger $l$ with increasing $k$, which is
not expected for a localized model. This is likely explained by the looseness
of the bound Eq. (B.2); while it is correct for the leading term in $D$, the
actual prefactors for terms of decreasing order in $D$ become increasingly
smaller than those in the given bound. We expect that a formula replacing Eq.
(B.2) in which the enumerative calculations are exact and the upper bound only
comes from $D$ would not possess this unexpected behaviour, although we have
been unable to derive such a formula. We nevertheless see the expected result
that with $k$ held constant we have an exponential decay of the ratio of
support for increasing $l$, and that increasing $D$ corresponds to an
increasing rate of decay.

## References

  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Edwards and Thouless [1972] J. T. Edwards and D. J. Thouless, Numerical studies of localization in disordered systems, J. Phys. C 5, 807 (1972).
  * Abrahams  _et al._ [1979] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Scaling theory of localization: Absence of quantum diffusion in two dimensions, Phys. Rev. Lett. 42, 673 (1979).
  * Abrahams [2010] E. Abrahams, ed.,  _50 Years of Anderson Localization_ (World Scientific, Singapore, 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 (2006).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. Prosen, and P. Prelovšek, Many-body localization in the heisenberg $xxz$ magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Pal and Huse [2010] A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field heisenberg chain, Phys. Rev. B 91, 081103 (2015).
  * Luitz  _et al._ [2016] D. J. Luitz, N. Laflorencie, and F. Alet, Extended slow dynamical regime close to the many-body localization transition, Phys. Rev. B 93, 060201 (2016).
  * Serbyn  _et al._ [2013] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013).
  * Vosk  _et al._ [2015] R. Vosk, D. A. Huse, and E. Altman, Theory of the many-body localization transition in one-dimensional systems, Phys. Rev. X 5, 031032 (2015).
  * Altman and Vosk [2015] E. Altman and R. Vosk, Universal dynamics and renormalization in many-body-localized systems, Annual Review of Condensed Matter Physics 6, 383 (2015).
  * Doggen  _et al._ [2018] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Many-body localization and delocalization in large quantum chains, Phys. Rev. B 98, 174202 (2018).
  * Sierant  _et al._ [2020] P. Sierant, D. Delande, and J. Zakrzewski, Thouless time analysis of anderson and many-body localization transitions, Phys. Rev. Lett. 124, 186601 (2020).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Potter  _et al._ [2015] A. C. Potter, R. Vasseur, and S. A. Parameswaran, Universal properties of many-body delocalization transitions, Phys. Rev. X 5, 031033 (2015).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).
  * Suntajs  _et al._ [2020a] J. Suntajs, J. Bonca, T. Prosen, and L. Vidmar, Quantum chaos challenges many-body localization, Phys. Rev. E 102, 062144 (2020a).
  * Suntajs  _et al._ [2020b] J. Suntajs, J. Bonca, T. Prosen, and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020b).
  * Kiefer-Emmanouilidis  _et al._ [2020a] M. Kiefer-Emmanouilidis, R. Unanyan, J. Sirker, and M. Fleischhauer, Bounds on the entanglement entropy by the number entropy in non-interacting fermionic systems, SciPost Phys. 8, 083 (2020a).
  * Kiefer-Emmanouilidis  _et al._ [2020b] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Evidence for unbounded growth of the number entropy in many-body localized phases, Phys. Rev. Lett. 124, 243601 (2020b).
  * Kiefer-Emmanouilidis  _et al._ [2021a] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Slow delocalization of particles in many-body localized phases, Phys. Rev. B 103, 024203 (2021a).
  * Kiefer-Emmanouilidis  _et al._ [2021b] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Unlimited growth of particle fluctuations in many-body localized phases, Annals of Physics , 168481 (2021b).
  * Kiefer-Emmanouilidis  _et al._ [2022a] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Particle fluctuations and the failure of simple effective models for many-body localized phases, SciPost Phys. 12, 034 (2022a).
  * Kiefer-Emmanouilidis  _et al._ [2022b] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Comment on "Resonance-induced growth of number entropy in strongly disordered systems", arXiv: 2203.06689 (2022b).
  * Sels and Polkovnikov [2021] D. Sels and A. Polkovnikov, Dynamical obstruction to localization in a disordered spin chain, Phys. Rev. E 104, 054105 (2021).
  * LeBlond  _et al._ [2021] T. LeBlond, D. Sels, A. Polkovnikov, and M. Rigol, Universality in the onset of quantum chaos in many-body systems, Phys. Rev. B 104, L201117 (2021).
  * Avdoshkin and Dymarsky [2020] A. Avdoshkin and A. Dymarsky, Euclidean operator growth and quantum chaos, Phys. Rev. Res. 2, 043234 (2020).
  * 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).
  * Heveling  _et al._ [2022] R. Heveling, J. Wang, and J. Gemmer, Numerically probing the universal operator growth hypothesis, Phys. Rev. E 106, 014152 (2022).
  * Fieker  _et al._ [2017] C. Fieker, W. Hart, T. Hofmann, and F. Johansson, Nemo/Hecke: Computer algebra and number theory packages for the Julia programming language, in  _Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation_, ISSAC ’17 (Association for Computing Machinery, New York, NY, USA, 2017) p. 157-164.
  * Dominici [2009] D. Dominici, Asymptotic analysis of the bell polynomials by the ray method, Journal of Computational and Applied Mathematics 233, 708 (2009), 9th OPSFA Conference.
  * Sels [2022] D. Sels, Bath-induced delocalization in interacting disordered spin chains, Physical Review B 106, L020202 (2022).
  * Aubry and André [1980] S. Aubry and G. André, Analyticity breaking and anderson localization in incommensurate lattices, Ann. Israel Phys. Soc. 3, 18 (1980).
  * Jitomirskaya [1999] S. Y. Jitomirskaya, Metal-insulator transition for the almost mathieu operator, Annals of Mathematics 150, 1159 (1999).
  * Bethe [1931] H. Bethe, Z. Phys. 71, 205 (1931).
  * Essler  _et al._ [2005] F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The One-Dimensional Hubbard model_ (Cambridge University Press, Cambridge, 2005).
  * Imbrie [2016a] J. Z. Imbrie, Diagonalization and many-body localization for a disordered quantum spin chain, Phys. Rev. Lett. 117, 027201 (2016a).
  * Imbrie [2016b] J. Imbrie, On Many-Body Localization for Quantum Spin Chains, Jour. Stat. Phys. 163, 998 (2016b).
  * Morningstar  _et al._ [2020] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Many-body localization near the critical point, Phys. Rev. B 102, 125134 (2020).
  * Imbrie  _et al._ [2017] J. Z. Imbrie, V. Ros, and A. Scardicchio, Local integrals of motion in many-body localized systems, Annalen der Physik 529, 1600278 (2017).
  * Lu  _et al._ [2023] B. Lu, C. Bertoni, S. J. Thomson, and J. Eisert, Measuring out quasi-local integrals of motion from entanglement (2023), arXiv:2301.01787 [cond-mat.dis-nn] .
  * Thomson and Schirò [2023] S. J. Thomson and M. Schirò, Local integrals of motion in quasiperiodic many-body localized systems, SciPost Phys. 14, 125 (2023).
  * Schrieffer and Wolff [1966] J. R. Schrieffer and P. A. Wolff, Relation between the anderson and kondo hamiltonians, Phys. Rev. 149, 491 (1966).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer-wolff transformation for quantum many-body systems, Annals of Physics 326, 2793 (2011).
  * Cao [2021] X. Cao, A statistical mechanism for operator growth, Journal of Physics A: Mathematical and Theoretical 54, 144001 (2021).
  * Sirker [2020] J. Sirker, Transport in one-dimensional integrable quantum systems, SciPost Phys. Lect. Notes , 17 (2020).
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).
  * Ha  _et al._ [2023] H. Ha, A. Morningstar, and D. A. Huse, Many-body resonances in the avalanche instability of many-body localization, Phys. Rev. Lett. 130, 250405 (2023).
  * Sels and Polkovnikov [2023] D. Sels and A. Polkovnikov, Thermalization of dilute impurities in one-dimensional spin chains, Phys. Rev. X 13, 011041 (2023).
  * Peacock and Sels [2023] J. C. Peacock and D. Sels, Many-body delocalization from embedded thermal inclusion, Physical Review B 108, 1020201 (2023).