# Probing localization properties of many-body Hamiltonians via an imaginary
vector potential

Liam O’Brien<EMAIL_ADDRESS>Department of Physics and Institute for
Quantum Information and Matter, California Institute of Technology, Pasadena,
CA 91125  Gil Refael<EMAIL_ADDRESS>Department of Physics and Institute
for Quantum Information and Matter, California Institute of Technology,
Pasadena, CA 91125

###### Abstract

Identifying and measuring the “localization length” in many-body systems in
the vicinity of a many-body localization transition is difficult. Following
Hatano and Nelson, a recent work Heußen et al. (2021) introduced an “imaginary
vector potential” to a disordered ring of interacting fermions, in order to
define a many-body localization length (corresponding, in the non-interacting
case, to the end-to-end Green’s function of the hermitian system). We extend
these results, by connecting this localization length to the length scale
appearing in the avalanche model of delocalization. We use this connection to
derive the distribution of the localization length at the MBL transition,
finding good agreement with our numerical observations. Our results
demonstrate how a localization length defined as such probes the localization
of the underlying ring, without the need to explicitly construct the l-bits.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Conventional statistical mechanics assumes that isolated systems reach thermal
equilibrium - effectively acting as their own heat bath - when left to evolve
under their own dynamics. The process by which this happens for quantum
systems is highly non-trivial Deutsch (1991); Srednicki (1994); Rigol et al.
(2008); Rigol and Srednicki (2012), and has been the subject of intense
research in recent decades - see Mori et al. (2018); Ueda (2020); Nandkishore
and Huse (2015) for recent reviews. Through this work, certain classes of
systems have been discovered which appear to violate this assumption, and do
not thermalize under their own dynamics. One such class that has received much
theoretical attention are lattice systems with quenched disorder. For non-
interacting systems, Anderson localization occurrs, and the single-particle
eigenstates are exponentially localized around a specific lattice site, with a
decay length that is a function of the disorder and the energy Anderson
(1958). Interacting systems display many-body localization (MBL), which is
characterized by zero DC transport, area-law entanglement entropy, Poissonian
level statistics, logarithmic growth of entanglement entropy, long-time memory
of initial conditions, and more Nandkishore and Huse (2015). MBL has emerged
as a platform for novel quantum order Huse et al. (2013); Pekker et al.
(2014); Chandran et al. (2014); Vasseur et al. (2015); Potter and Vishwanath
(2015); Bahri et al. (2015); Vasseur et al. (2016); Kells et al. (2018);
Iadecola and Schecter (2018); Friedman et al. (2018); Kuno (2019) (see
Parameswaran and Vasseur (2018) for a review), with implications for quantum
information Santos et al. (2005); Altshuler et al. (2010); Laumann et al.
(2015); Yao et al. (2015); Khemani et al. (2015); Goihl et al. (2020).

The existence of a many-body localized regime has been demonstrated
numerically in small systems Pal and Huse (2010); Bauer and Nayak (2013);
Žnidarič et al. (2008); Kjäll et al. (2014); Luitz et al. (2015) and proven
analytically Imbrie (2016) for one-dimensional spin chains under certain
assumptions, but many open questions still remain. In particular, the nature
(universality class) of the transition between the thermal/ergodic phase and
the MBL phase still remains to be understood Nandkishore and Huse (2015). Some
recent works even call into question whether the transition exists at all in
the thermodynamic limit Šuntajs et al. (2020); Sels and Polkovnikov (2021);
subsequent works have argued the observed effects underlying those conclusions
are the product of finite-size effects Abanin et al. (2021) and many-body
resonances Crowley and Chandran (2022); Morningstar et al. (2022). Much
remains to be done in probing the nature and stability of the localization in
systems that are believed to be MBL.

The most natural way to quantify localization is through the localization
length. In the non-interacting case, this is straightforward: it is the decay
length of the single-particle eigenstates. In an interacting system, the
desired “localization length” is less obvious, since there is no notion of
single particle orbitals. One of the most useful ways to quantify the
localization, then, is through the use of localized conserved quantities. In
addition to the properties mentioned above, the MBL phase comes equipped with
a complete set of (quasi-)local integrals of motions (LIOMs) Huse et al.
(2014); Ros et al. (2015) \- the so-called “l-bits”. Roughly speaking, the
l-bits are conserved quantities that are obtained by “dressing” the physical
degrees of freedom (“p-bits”) with unitaries that are exponentially localized
in real space. The decay of the l-bits thus offers a way of calculating a
localization length in interacting systems.

A number of methods have been put forth to address the issue of constructing
the l-bits Serbyn et al. (2013); Chandran et al. (2015); Rademaker and Ortuño
(2016); Pekker et al. (2017), but the procedure is computationally taxing and
physically ambiguous (since the mapping of the set of l-bit configurations
onto the computational basis states is not unique). This ambiguity means that
the localization lengths measured depend on how the l-bits are constructed. In
principle, it is possible to uniquely specify an assignment of l-bits that is
“most localized” in the original basis Huse et al. (2014), but there is no
algorithm known rigorously to construct such an assignment (though algorithms
based on Wegner-Wilson flow Wegner (1994); Kehrein (2006); Pekker et al.
(2017); Quito et al. (2016) work fairly well in practice). In spite of these
difficulties, several numerical studies Peng et al. (2019); Kulshreshtha et
al. (2018); Varma et al. (2019) have extracted localization lengths from the
l-bit decay and analyzed properties of the finite-size MBL phase and MBL-
thermal crossover.

An alternate method for finding a localization length is to introduce an
imaginary vector potential, or tilt, that makes hopping preferential in one
direction Heußen et al. (2021). By varying the strength of the vector
potential/tilt, one finds that various eigenvalues will develop non-zero
imaginary parts (as the system is no longer Hermitian). The point at which a
given eigenvalue develops a non-zero imaginary part in this modified model
defines a length scale $\xi$ that, in non-interacting systems, directly
measures the localization length (defined in terms of an end-to-end Green’s
function) of the analogous (single-particle) eigenstate of the underlying
model without the vector potential. This procedure does not require
constructing the l-bits, and Ref Heußen et al. (2021) showed the corresponding
$\xi$ can be used to identify the critical parameters of the MBL-thermal
crossover in interacting systems.

Here, we connect the length scale $\xi$ to the avalanche model of
delocalization. Much of the current literature describes the asymptotic MBL to
thermal transition in terms of such “avalanches”, whereby rare (locally)
thermal regions are able to induce a cascade of thermalization that overpowers
the localization in the rest of the system De Roeck and Huveneers (2017);
Luitz et al. (2017); Thiery et al. (2018, 2017); Johri et al. (2015). The
instability of localized systems to avalanches is characterized by a length
scale $\lambda$, which captures the decay of matrix elements coupling the
l-bits and rare thermal regions. It is natural, then, to expect that $\xi$ and
$\lambda$ are related to each other, as they both act as a measure for the
localization-delocalization transition; we make this relation explicit.

Moreover, $\xi$ has a distribution (with respect to quenched disorder). We
calculate this distribution numerically at the finite-size MBL-thermal
crossover, and examine whether the relation between $\xi$ and $\lambda$ holds
for the distributions as well. Indeed, we find it does, so that the
distributions of $\xi$ at the MBL crossover contain information about the
decay of matrix elements at the corresponding transition.

Our paper is organized as follows: In section II, we present and discuss the
model we will study and some of its salient features. In section III we show
$\xi$ can be predicted by proliferation of non-hermitian avalanches, and we
connect $\xi$ to the decay length $\lambda$ appearing in the avalanche picture
of the delocalization transition. Finally, in section IV, we use this
connection to derive an analytic form of the distribution for $\xi$, and
compare with the numerically observed histograms.

##  II Non-hermitian MBL Model

###  II.1 Hamiltonian

We consider the following Hamiltonian for spinless fermions on a one-
dimensional lattice with sites $i=1,\ldots,L$:

$\displaystyle
H=t\sum_{i=1}^{L}[e^{g}c_{i}^{{\dagger}}c_{i+1}+e^{-g}c_{i+1}^{{\dagger}}c_{i}]+\sum_{i=1}^{L}h_{i}n_{i}+U\sum_{i=1}^{L}n_{i}n_{i+1}$

(1)

where $c_{i}^{{\dagger}}/c_{i}$ are fermionic creation/annhiliation operators,
$n_{i}=c_{i}^{{\dagger}}c_{i}$ number operators, $t\in\mathbb{R}$ is the
hopping amplitude, $g\geq 0$ is the “imaginary vector potential” or “tilt”,
$U>0$ is a parameter denoting the strength of interactions, and the $h_{i}$’s
are random variables drawn independently and identically from a distribution
characterized by a “disorder strength” $W$. This Hamiltonian describes
particles hopping on a disordered lattice, subject to nearest neighbor
repulsion, with an imaginary vector potential (of magnitude $g$) making left
hopping preferential - see Fig. 1 for a schematic of this model.

Figure 1: A schematic of the model Hamiltonian (1). Fermionic particles
(filled red circles) can hop (preferentially to the left) with strength
$te^{\pm g}$, interact with their nearest neighbor with strength $U$, and are
also subject to a random potential whose strength is determined by the
disorder width $W$.

The Hamiltonian in (1) has been studied as an effective model to describe
driven open systems Panda and Banerjee (2020), and can also be mapped onto the
statistical mechanics of depinning flux lines from columnar defects in two-
dimensional type II superconductors via the path integral formalism Hatano and
Nelson (1996); Hatano and Nelson (1998); Hatano and Nelson (1997). In such a
mapping, $g$ characterizes an external magnetic field that is tilted with
respect to the random columnar defects - hence the choice of names for $g$.

Since $g$ is a real number, we generically expect some eigenvalues to be
complex for $g\neq 0$. The values of $g$ at which certain eigenvalues develop
a non-zero imaginary part are examples of so-called exceptional points Kato
(1995), which have been a topic of research since the 1990s Heiss and Sannino
(1990); Heiss and Sannino (1991); Heiss (2000, 2001); Heiss and Harney (2001);
Heiss (2004, 2008, 2012); Luitz and Piazza (2019). Of particular interest has
been the connection between the exceptional points and
localization/delocalization of eigenstates in models similar to (1) Hatano and
Nelson (1996); Hatano and Nelson (1998); Hatano and Nelson (1997); Hamazaki et
al. (2019); Feinberg and Zee (1999a); Feinberg and Zee (1999b); Brézin and Zee
(1998). We seek to exploit this connection, as we lay out in the next section.

###  II.2 Exceptional points of $H$ and their connection to
localization/delocalization

To expand on the connection between exceptional points and delocalization in
the model (1), it is instructive to consider the following “gauge
transformation” Heußen et al. (2021):

$S=\exp\left(\sum_{j=1}^{L}jgn_{j}\right)$ (2)

Applying this transformation to our Hamiltonian (1) with open boundary
conditions ($c_{L+1}=0$) eliminates $g$ entirely. The Hamiltonian is similar
to a Hermitian operator, and thus has a real spectrum, for any $g$. Put
another way, we can “gauge away” the imaginary vector potential. Conversely,
if we have periodic boundary conditions ($c_{L+1}=c_{1}$) the presence of $g$
is not removed, but rather is shifted entirely to the bond between sites $L$
and $1$ (or any other two sites, by an appropriate shift of indices in (2)).
In this case, the non-hermiticity cannot be removed; there is a fixed “flux”
$iLg$ through the ring. See Fig. 2(a) for an illustration.

(a)

(b)

Figure 2:  a) Illustration of the imaginary flux $\Phi=iLg$ bound through our
system (with periodic boundary conditions), as well as the action of the gauge
transformation (2). b) Illustration of the competition between the decay of
the wavefunction $\left\lvert\psi_{\epsilon}\right\rangle$ and the imaginary
flux $iLg$. Heuristically, we expect the exceptional point $g_{c}(\epsilon)$
to occur where the product of the end-to-end tunneling amplitude $\left\langle
1\right\rvert\left.\psi_{\epsilon}\right\rangle\left\langle\psi_{\epsilon}\right\rvert\left.L\right\rangle$
and hopping strength $|t|e^{Lg}$ is comparable to the level spacing.

Now, let us consider the single-particle sector of the Hamiltonian (1) with
open boundary conditions and $g=0$; we know that the states will be localized
for any non-zero disorder strength $W$. Let us consider an eigenstate
$\left\lvert\psi_{\epsilon}\right\rangle$ at energy density $\epsilon$ 111For
a given eigenvalue $E$, the energy density is defined as the excitation energy
of $E$ normalized by the bandwidth:
$\epsilon=\frac{E-E_{\min}}{E_{\max}-E_{\min}}$ Here $E_{\max}$, $E_{\min}$
are maximum and minimum energies in the spectrum, respectively. in this open
chain. This eigenstate will be localized at some lattice site - call it $j$.
Based on the discussion of gauge freedom above, let us introduce an imaginary
flux by adding a term of the form
$e^{Lg}c_{L}^{{\dagger}}c_{1}+e^{-Lg}c_{1}^{{\dagger}}c_{L}$ to our system.
Borrowing intuition from perturbation theory, we expect the non-hermiticity to
become important when:

$\left\langle\psi_{\epsilon}\right\rvert|t|e^{Lg}c_{L}^{{\dagger}}c_{1}\left\lvert\psi_{\epsilon}\right\rangle\sim\delta$
(3)

where $\delta$ is the level spacing in the chain. See Fig. 2(b) for an
illustration.

Using the localized state ansatz
$\left\lvert\psi_{\epsilon}\right\rangle\sim\sum_{k}e^{-|k-j|/\xi_{\epsilon}}\left\lvert
k\right\rangle$ in (3) above, one finds the exceptional point
$g_{c}(\epsilon)$ for this eigenstate to be $g_{c}(\epsilon)\sim
1/\xi_{\epsilon}$ to leading order in $L$. In other words, the tilt competes
directly with the localization of the underlying hermitian Hamiltonian. This
fact has been studied in detail by Hatano and Nelson for the non-interacting
case Hatano and Nelson (1996); Hatano and Nelson (1998) and Hamazaki et al
Hamazaki et al. (2019) for the interacting case. Hamzaki et al focused on the
localization properties of the non-hermitian Hamiltonian (1) at nonzero $g$,
as opposed to using the real-complex transition to probe the localization
properties of the underlying $g=0$ Hamiltonian (which is what the relation
$g_{c}\sim 1/\xi$ encodes).

The above argument, while suggestive, can be made precise. Specifically, let
$\left\lvert\psi_{n}\right\rangle$ denote the $n^{\text{th}}$ single-particle
eigenstate of the $g=0$ Hamiltonian (1) with open boundary conditions, and let
$g_{c,\,n}$ be the exceptional point for the $n^{\text{th}}$ eigenstate in the
closed chain with tilt (ordered by real part of the eigenvalue). Then one can
show Heußen et al. (2021):

$g_{c,\,n}=\frac{1}{\xi_{n}}+\mathcal{O}\left(\frac{1}{L}\right)$ (4)

where $1/\xi_{n}=-\frac{1}{L}\log(\left\langle\psi_{n}\right\rvert
c_{1}^{{\dagger}}c_{L}\left\lvert\psi_{n}\right\rangle)$ is the “eigenstate
localization length” (equivalently, the logarithm of the end-to-end Green’s
function in the vicinity of $E_{n}$). While derived only for the single-
particle sector of (1), Ref Heußen et al. (2021) showed that defining a length
scale via (4) in interacting systems has the expected behavior of a
localization length - namely, a finite-size scaling collapse of $\xi$
identifies critical parameters that agree with other numerical studies of
finite-size MBL systems, and $\xi$ appears to diverge at the asymptotic
transition. Thus, by searching for the exceptional points of the Hamiltonian
(1) as we vary $g$, we can directly probe the localization properties of the
underlying $g=0$ Hamiltonian.

##  III Connection to Avalanche Model

In this section, we seek to connect the localization length $\xi$ to the
length scale $\lambda$ appearing in the avalanche model of delocalization. We
do so in several steps.

###  III.1 Hermitian Avalanche Model

We begin with a brief review of the avalanche model of the localization-
delocalization transition De Roeck and Huveneers (2017); Luitz et al. (2017);
Thiery et al. (2018, 2017); Johri et al. (2015). In this picture, localized
chains are coupled weakly to a thermal bath (or thermal sub-region of the
chain), which can thermalize the sites bordering it. The new, larger,
“effective” thermal bath comprised of the original bath and its thermalized
neighbors can now thermalize the subsequent sites bordering it, and this can
continue, triggering an “avalanche” that will either halt at some point -
leaving the chain (partially) localized - or thermalize the whole chain.

To be more concrete, we consider the following toy model that captures the
essential physics of such avalanches De Roeck and Huveneers (2017); Luitz et
al. (2017). The system is a 1d chain of length $L+L_{b}$, with sites
$-L_{b}+1,\ldots,0$ being a thermal “bath”, and sites $1,\ldots,L$ being many-
body localized l-bits. The Hamiltonian is

$H=H_{\text{bath}}+\mathcal{G}_{0}\sum_{j=1}^{L}e^{-j/\lambda}[b_{0}^{{\dagger}}\tilde{c}_{j}+\tilde{c}_{j}^{{\dagger}}b_{0}]+\sum_{j=1}^{L}u_{j}\tilde{n}_{j}$
(5)

where $H_{\text{bath}}$ is the bath Hamiltonian (acting only on the bath
sites), $b_{i}/b_{i}^{\dagger}$ are bath creation/annihilation operators,
$\tilde{c}_{i}/\tilde{c}_{i}^{\dagger}/\tilde{n}_{i}$ are l-bit
creation/annihilation/number operators, and $\\{u_{i}\\}_{i=1}^{L}$ are the
single-particle energies associated with each l-bit. Additionally,
$e^{-1/\lambda}\in(0,1]$ is the base of exponential decay of the matrix
elements coupling the bath to the localized part of the chain, with $\lambda$
being the decay length of said matrix elements. This model neglects
interactions between the l-bits, which are argued to induce only higher order
corrections Luitz et al. (2017).

In the avalanche picture of delocalization, the eigenstate thermalization
hypothesis (ETH) holds for a given l-bit if and only if it sucessfuly
hybridizes with the bath according to Fermi’s golden rule Potirniche et al.
(2019). Symbolically, this criteria amounts to the condition
$\mathcal{T}\gg\delta$, where $\delta$ is the level spacing in the bath, and
$\mathcal{T}$ is a matrix element of the hopping
$b_{0}^{\dagger}\tilde{c}_{i}$ between eigenstates of the unperturbed
Hamiltonian $H_{\text{bath}}+\sum_{i=1}^{L}u_{i}\tilde{n}_{i}$. Assuming ETH
for the bath, a typical matrix element between two energy eigenstates has the
form $\mathcal{T}\sim\kappa\sqrt{\rho(\epsilon,\omega)\delta}$, where the
spectral function $\rho$ is a smooth positive function, $\epsilon$ and
$\omega$ are the average and difference of the bath energies of the
eigenstates, and $\kappa=\mathcal{G}_{0}e^{-j/\lambda}$ is direct coupling
between the bath and l-bit $j$ De Roeck and Huveneers (2017).

The l-bit most strongly coupled to the bath will be at site 1, and this l-bit
hybridizes with the bath if
$\mathcal{G}_{0}e^{-1/\lambda}\sqrt{\rho(\epsilon,\omega)/\delta}\gg 1$.
Should this l-bit hybridize, it will effectively become part of the bath; this
(approximately) halves the level spacing $\delta\mapsto\delta/2$ and
(approximately) leaves the spectral function unchanged:
$\rho^{\prime}(\epsilon,\omega)\approx\rho(\epsilon,\omega)$. One can now pose
the hybridization criteria to the l-bit at site 2, with the level spacing of
the “effective” bath half that of the original bath level spacing. Should the
l-bit at site 2 hybridize with the bath, the level spacing will be further
reduced by a factor of 2. This process continues at each site a distance $r$
from the bath so long as the hybridization condition

$\frac{\mathcal{T}(r)}{\delta(r)}\gg
1\quad\Longrightarrow\quad\mathcal{G}_{0}e^{-r/\lambda}2^{r/2}\sqrt{\frac{\rho(\epsilon,\omega)}{2\delta_{0}}}\gg
1$ (6)

holds. Here, we’ve taken $\delta(r)\approx\delta_{0}2^{-r+1}$, with
$\delta_{0}$ the level spacing of the original bath.

For the entire chain to thermalize, (6) should hold for all $r$. Since the
left hand side of (6) is monotonic, full thermalization can be determined by
the behavior at the end of the chain ($r=L$). Taking logarithms, full
thermalization of the chain thus amounts to:

$-\frac{1}{\lambda}+\frac{\log
2}{2}+\mathcal{O}\left(\frac{1}{L}\right)\gtrapprox 0$ (7)

In the thermodynamic limit, we can see this condition will always be satisfied
for $\lambda>2/\log 2$. We thus see the length scale $\lambda$ controls the
delocalization transition, which occurs at a critical decay length
$\lambda_{c}=2/\log 2$.

###  III.2 Non-hermitian delocalization & avalanches

In numerical studies of non-hermitian Hamiltonians with tilt $g$ whose $g=0$
counterpart is many-body localized, it has been observed that a non-hermitian
MBL regime exists for small $g$, and is eventually destroyed for large enough
$g$ Hamazaki et al. (2019); Panda and Banerjee (2020). We conjecture that this
non-hermitian delocalization occurrs by the same avalanche mechanism as in the
hermitian case. That is, we propose that the imaginary vector potential helps
couple ergodic grains to the rest of the chain in such a way as to
successively delocalize neighboring sites of the chain.

Note that Ref Hamazaki et al. (2019) also argues that the real-complex
transition roughly coincides with the non-hermitian delocalization transition
in their numerics, and that the coincidence is exact in the thermodynamic
limit. Under our conjecture, the exceptional points $g_{c}$ thus measure the
location of this avalanche-based transition (i.e.,
$g_{c}\sim\lambda_{c}-\lambda$).

The conjecture that a non-hermitian avalanche drives the non-hermitian
delocalization transition is a highly non-trivial statement, so we first
verify it on a simple toy model. To that end, we consider the following non-
hermitian version of the avalanche Hamiltonian:

$H=H_{\text{bath}}+\sum_{j=1}^{L}u_{j}\tilde{n}_{j}+\mathcal{G}_{0}\sum_{j=1}^{L}e^{-j/\lambda}\left[e^{jg}b_{0}^{{\dagger}}\tilde{c}_{j}+e^{-jg}\tilde{c}_{j}^{{\dagger}}b_{0}\right]+\,\mathcal{G}_{0}\sum_{j=1}^{L}e^{-j/\lambda}\left[e^{jg}\tilde{c}_{L-j+1}^{{\dagger}}b_{-L_{b}+1}+e^{-jg}b_{-L_{b}+1}^{{\dagger}}\tilde{c}_{L-j+1}\right]$
(8)

This is the Hamiltonian (5), with the inclusion of a non-hermitian tilt. Note
also that we include hopping from both ends of the bath, since the effect of
$g$ can be gauged away for open boundary conditions.

Figure 3 shows numerical results from diagonalizing the Hamiltonian (8) for
several system sizes, in the half-filling sector, for $e^{-1/\lambda}=0.3$. We
can see a crossover in both the half-chain entanglement entropy and the
inverse participation ratio from constant to increasing/decaying with system
size (respectively), indicating a localization to delocalization crossover.
These crossovers also appear to occurr at the value of $g$ we would expect
based on modifying the avalanche criteria (6) with $e^{-r/\lambda}\mapsto
e^{-r/\lambda+rg}$, i.e. $g=1/\lambda-1/\lambda_{c}$ (with $\lambda_{c}=2/\log
2$). We see similar crossovers, also at the expected values of $g$, if we vary
the value of $e^{-1/\lambda}$ \- see appendix C for additional examples. These
results support the idea that a non-hermitian avalanche mechanism is driving
delocalization in this model.

Figure 3(c) also shows the fraction $f_{\text{comp}}$ of disorder realizations
that have a non-real eigenvalue at energy density $\epsilon=0.5$. Though not
exactly coinciding with the expected transition value of
$g=1/\lambda-1/\lambda_{c}$, the window over which a non-zero fraction of
disorder realizations have a complex eigenvalue overlaps with the crossover
window observed in Figs. 3(a) and 3(b). This suggests that the real-complex
transition may not only coincide with the many-body localization-
delocalization transition (as previously suggested in the literature), but
that it is consistent with a non-hermitian avalanche mechanism.

(a)

(b)

(c)

Figure 3: Disorder-averaged localization metrics of the eigenstates of the
non-hermitian avalanche Hamiltonian (8) for $e^{-1/\lambda}=0.3$ and various
$L$, as a function of $g$. Going clockwise: a) half-chain entanglement entropy
$\langle S_{E}\rangle$, b) Inverse participation ratio $\langle IPR\rangle$,
and c) Fraction $f_{\text{comp}}$ of disorder realizations whose central
eigenvalue ($\epsilon=0.5$) is non-real (note that
$\left\langle\cdot\right\rangle$ denotes a disorder average). From a) and b),
we see a crossover from localized behavior to delocalized behavior, and from
c), we see a similar real-complex crossover. Insets show a “zoomed in” look at
the crossover region. These crossovers appear to occurr at
$g=1/\lambda-1/\lambda_{c}$ (black dashed lines), the value one would expect
from a non-hermitian avalanche criterion, suggesting non-hermitian avalanches
are responsible for these crossovers. All quantities were computed for the
eigenstate at the center of the spectrum ($\epsilon=0.5$) in the half-filling
sector with $\mathcal{G}_{0}=1$, $L_{b}=3$, and $H_{\text{bath}}$ an $8\times
8$ GOE matrix. The l-bit energies were drawn from a uniform distribution
$u_{i}\sim\text{Uni}[-w,w]$ with $w=10$, and we averaged all quantities over
1000 ($L=5,7$) and 500 ($L=9$) disorder realizations. The entanglement
entropies were computed using right eigenvectors only.

###  III.3 Mapping chain hamiltonians to avalanche hamiltonians

Having evidence that a non-hermitian avalanche drives the delocalization and
real-complex transition in a simple toy model, we now would like to verify
that the exceptional points $g_{c}(\epsilon)$ in more generic Hamiltonians
such as (1) are described by a similar mechanism. To do so, we need a way to
explicitly connect such Hamiltonians to “avalanche” Hamiltonians resembling
(8).

Let us consider a generic system comprised of a thermal bath/grain and a chain
(which we take to be MBL). The Hamitonian is of the form:

$H=H_{\text{bath}}+H_{\text{chain}}+H_{\text{bc}}$ (9)

where $H_{\text{bath}}$ and $H_{\text{chain}}$ act only on the bath and chain
(respectively), and $H_{\text{bc}}$ couples the bath and chain subsystems. To
bring this into the desired “avalanche” form, we need to re-express
$H_{\text{chain}}$ and $H_{\text{bc}}$ in terms of $l$-bits. In principle,
this can be done straightforwardly by diagonalizing $H_{\text{chain}}$ and
performing a change of basis on $H_{\text{bc}}$. When using a generic
diagonalization routine, however, it is not obvious how to extract the l-bit
occupation numbers for each eigenstate. As the avalanche model relies on a
cascade effect from thermalizing successive l-bits, we need to employ a
diagonalization routine that allows us to access this crucial information.

To achieve this, we employ a mapping that combines the “displacement
transformations” of Rademaker and Ortuño Rademaker and Ortuño (2016) with the
principles of Wegner-Wilson Flow Pekker et al. (2017); Quito et al. (2016).
Specifically, we construct and apply displacement transformations to eliminate
individual off-diagonal terms in the Hamiltonian
$H_{0}=H_{\text{bath}}+H_{\text{chain}}$. We iterate this procedure,
repeatedly eliminating the largest remaining off-diagonal term (as in the flow
equation approach) until all such terms in $H_{0}$ are below some tolerance -
see Appendix D for details. At the end, this procedure yields similarity
transformations $U_{\text{bath}}$ and $U_{\text{chain}}$ that diagonalize
$H_{\text{bath}}$ and $H_{\text{chain}}$, respectively. We then apply these
transformations to $H_{\text{bc}}$, which generates an avalanche-like series
of terms involving hopping between bath eigenstates and the l-bits. Figure 4
illustrates this procedure schematically.

(a)

(b)

Figure 4: A schematic illustration of how we map chains with a bath to
avalanche-like Hamiltonians. a) We first diagonalize the decoupled chain and
bath system (above the arrow), by iteratively applying displacement
transformations to eliminate the largest off-diagonal hopping terms (solid
lines). This leaves the system in diagonal form (below arrow), with
$\left\lvert\psi_{\text{b}}^{i}\right\rangle$ the bath eigenstates and squares
representing the l-bits $\ell_{i}$. b) We then apply the same transformations
found in (a) to the full Hamiltonian $H$, including the bonds between the bath
and chain (red lines, above arrow). This yields avalanche-like hopping between
the bath eigenstates and l-bits (red lines, below arrow).

There are several important points about this procedure that we wish to
emphasize. The first is that the transformation $U_{\text{chain}}$ obtained
from the Wegner-Wilson-flow like approach will well approximate the “near-
optimal” mapping of basis states to l-bit configurations alluded to in Section
I Pekker et al. (2017). Thus, we can read off the l-bit occupation numbers
from the physical site occupation numbers in the original computational basis.
The second is that by completely diagonalizing the decoupled bath/chain system
$H_{0}$ \- rather than diagonalizing only $H_{\text{chain}}$ \- we have access
to the spectrum of the decoupled system. We can therefore label each
l-bit/eigenstate hopping by its energy density in the uncoupled system,
allowing for direct comparison with $g_{c}$ at the same energy density. The
final key point is that, due to the gauge freedom, we can take $H_{0}$ to be
hermitian without loss of generality. More specifically, we can gauge all of
the flux onto the coupling Hamiltonian $H_{\text{bc}}$ (i.e., the red bonds in
Fig 4(b)), and perform our mapping on $H_{0}$ with $g=0$. The transformations
$U_{\text{bath}}$ and $U_{\text{chain}}$ will be unitary, and significantly
easier to construct. The influence of the imaginary flux is also now simply a
multiplicative weight for the various terms generated from $H_{\text{bc}}$ by
the change of basis (rather than being part of the transformations
$U_{\text{bath}}$ and $U_{\text{chain}}$, as it would be otherwise). We will
leverage these last two points in the next two sections.

###  III.4 Applying the mapping

Let us now apply the mapping described in the previous section to a concrete
model describing our system of interest. To that end, we take
$H_{\text{bath}}$ to be a $2^{L_{\text{b}}}\times 2^{L_{\text{b}}}$ GOE matrix
(where $L_{\text{b}}$ is the bath size) and

$\begin{split}H_{\text{chain}}=&\displaystyle\sum_{i=1}^{L-1}\left[c_{i}^{{\dagger}}c_{i+1}+c_{i+1}^{{\dagger}}c_{i}\right]\\\
&+\sum_{i=1}^{L}w_{i}n_{i}+U\sum_{i=1}^{L-1}n_{i}n_{i+1}\end{split}$ (10a)
$\begin{split}H_{\text{bc}}=&b_{0}^{{\dagger}}c_{1}+c_{1}^{{\dagger}}b_{0}+c_{L}^{{\dagger}}b_{-L_{b}+1}+b_{-L_{b}+1}^{{\dagger}}c_{L}\end{split}$
(10b)

where $b_{i}$, $c_{i}$ are bath/chain fermion operators (respectively),
$n_{i}$ are chain number operators, $U>0$ is the interaction strength, the
$w_{i}$ are independently and identically drawn from a distribution
characterized by strength $W$, and $L$ is the length of the chain. Note that
we have taken the bare hopping strength $t$ (c.f. the Hamiltonian (1)) to be
unity, so that all quantities are in units of $t$. Additionally, per the
discussion of the previous section, we apply the mapping to a hermitian
Hamiltonian without the vector potential $g$; we can restore the influence of
$g$ by placing a flux $i(L+1)g$ onto the bonds in and out of the bath.

For each disorder realization, after applying the mapping, we need to extract
the decay of matrix elements coupling each l-bit to the bath. To do so, we
first separate the bath-chain coupling $H_{\text{bc}}$ into couplings from the
left and right end of the bath:
$H_{\text{bc}}^{L}=c_{L}^{{\dagger}}b_{-L_{b}+1}+b_{-L_{b}+1}^{{\dagger}}c_{L}$
and $H_{\text{bc}}^{R}=b_{0}^{{\dagger}}c_{1}+c_{1}^{{\dagger}}b_{0}$,
respectively. In terms of the uncoupled eigenstates
$\left\lvert\Psi_{ij}\right\rangle=\left\lvert\phi_{i}\right\rangle\left\lvert\psi_{j}\right\rangle$
(where $\left\lvert\phi_{i}\right\rangle$, $\left\lvert\psi_{j}\right\rangle$
are eigenstates of $H_{\text{bath}}$ and $H_{\text{chain}}$, respectively), we
obtain matrix elements:

$\mathcal{A}^{ij;\,kl;\,s}=\left\langle\Psi_{ij}\right\rvert
H_{\text{bc}}^{s}\left\lvert\Psi_{kl}\right\rangle$

We wish to extract specific matrix elements for each l-bit $\ell$ and energy
density $\epsilon$. To that end, we choose $k=k_{\epsilon}$ and
$l=l_{\epsilon}$ above such that
$\left\lvert\Psi_{k_{\epsilon}l_{\epsilon}}\right\rangle$ is the eigenstate
closest to energy density $\epsilon$ in the uncoupled system. Then, for each
l-bit $\ell$, we define the amplitude $\mathcal{A}_{\ell}(\epsilon)$ via:

$\displaystyle\mathcal{A}_{\ell}(\epsilon)=\begin{cases}\max_{i}\left\lvert\left\langle\phi_{i}\right\rvert\left(\left\langle\psi_{l_{\epsilon}}\right\rvert\tilde{c}^{{\dagger}}_{\ell}\right)H^{R}_{\text{bc}}\left\lvert\Psi_{k_{\epsilon}l_{\epsilon}}\right\rangle\right\rvert&\ell\text{
occupied}\\\
\vspace{1pt}\max_{i}\bigg{\lvert}\left\langle\phi_{i}\right\rvert\left(\left\langle\psi_{l_{\epsilon}}\right\rvert\tilde{c}_{\ell}\right)H^{L}_{\text{bc}}\left\lvert\Psi_{k_{\epsilon}l_{\epsilon}}\right\rangle\bigg{\rvert}&\ell\text{
unoccupied}\end{cases}$

(11)

In words, we choose the largest (by amplitude) matrix element conecting
$\left\lvert\Psi_{k_{\epsilon}l_{\epsilon}}\right\rangle$ to a state with the
occupation of l-bit $\ell$ flipped and all other l-bit occupations the same.
Note that the choice of $H^{R}_{\text{bc}}$ versus $H^{L}_{\text{bc}}$
reflects the fact that hopping left (right) is enhanced (suppressed) when we
eventually restore the imaginary vector potential $g$; hence, we choose the
matrix element that will be enhanced by a factor of $e^{(L+1)g}$ in each case
above.

In the event that it is impossible to find such a hopping (if, for example,
l-bit $\ell$ is unoccupied but the bath is empty), then we examine the
eigenstate with energy density next closest to $\epsilon$.

In the context of an avalanche mechanism, we are interested in studying how
the couplings $\mathcal{A}_{\ell}(\epsilon)$ decay for each energy density. To
that end, we assume without loss of generality that the
$\mathcal{A}_{\ell}(\epsilon)$ are sorted in $\ell$ in descending order. The
amplitudes $\mathcal{A}_{\ell}(\epsilon)$, as well as their decay in $\ell$,
will generically be random (as they depend on the random disorder) - this is
in stark contrast with the deterministic decay of the canonical avalanche
model discussed in Section III.1. We do, however, observe that upon averaging
over disorder realizations, the amplitudes decay exponentially in $\ell$ \-
see Fig. 5 for examples of how this is borne out in chains of size $L=11$,
$L_{\text{b}}=3$, in the $N=3$ occupation sector 222The mapping described in
section III.3 is computationally very expensive to run, so we use a smaller
filling factor here in order to reach larger system sizes while keeping the
Hilbert space dimension roughly fixed. We have verified that we attain
qualitatively similar results at half-filling with smaller system sizes..
Consequently, the canonical avalanche model of Section III.1 emerges from our
results upon disorder averaging.

(a)

(b)

(c)

(d)

(e)

Figure 5:  Decay of $\left\langle\mathcal{A}_{\ell}(\epsilon)\right\rangle$
(with $\left\langle\cdot\right\rangle$ denoting disorder average) versus
$\ell$ at the center of the spectrum ($\epsilon=0.5$), for interaction
strengths a) $U=1$, b) $U=2$, c) $U=3$, d) $U=4$, e) $U=5$. Within each panel,
each line traces out the decay of
$\left\langle\mathcal{A}_{\ell}(\epsilon)\right\rangle$ for a given disorder
strength $W$, with the color of the line indicating the value of $W$. We see
that the average amplitudes appear to decay exponentially for all disorder
strengths, with the rate of decay increasing with the disorder strength. All
quantities are computed with $L=11$, $L_{b}=3$, in the $N=3$ sector, and the
averages are over 500 disorder realizations.

This mapping allows us to go beyond studying the disorder-averaged amplitudes
- in particular, we have access to the distributions of
$\mathcal{A}_{\ell}(\epsilon)$ and their decay in $\ell$. To that end, we
define the quantities

$\chi_{\ell}(\epsilon)=\frac{\mathcal{A}_{\ell}(\epsilon)}{\mathcal{A}_{\ell-1}(\epsilon)}$
(12)

, where we take $\mathcal{A}_{0}(\epsilon)=1$. Clearly then,
$\mathcal{A}_{\ell}(\epsilon)=\prod_{i=1}^{\ell}\chi_{i}(\epsilon)$, with
$0\leq\chi_{i}(\epsilon)\leq 1$ for each $\ell$ (as the
$\mathcal{A}_{\ell}(\epsilon)$ are sorted to be monotonically decreasing in
$\ell$). Such a decomposition amounts to replacing the deterministic decay
$e^{-1/\lambda}$ in (8) by the random, energy dependent
$\chi_{\ell}(\epsilon)$. The distributions of these $\chi_{\ell}(\epsilon)$
show very interesting behavior as we bring the system from MBL to ergodic.
Fig. 6 shows the observed probability densities of $\chi_{\ell}(\epsilon)$
(with $\ell=4$) as we tune the disorder strength from deep in the MBL regime
to ergodic. In the MBL and thermal regimes, the distributions are peaked near
0 and 1, respectively. This is consistent with the intuition that the
amplitudes should decay slowly in the thermal regime and quickly in the MBL
regime (which we observe to be the case on average in Fig. 5). Significantly,
there is a window between the two extreme regimes where the distributions are
approximately uniform. To make this statement more precise, we “fit” the
numerically observed probability distributions in the vicinity of the MBL-
thermal crossover, by minimizing the Akaike information over a large class of
known continuous probability distributions. The Akaike information is defined
as:

$I_{\text{A}}=2k-\log\hat{L}$ (13)

where $k$ is the number of parameters in the model, and $\hat{L}$ is the
likelihood function. We choose to minimize the Akaike information, as it acts
as an unbiased estimator for the expected Kullback-Leibler divergence (whereas
the Bayesian information is a biased estimator), a standard measure for the
information difference between two probability density functions Burnham and
Anderson (2004). Fig 7 shows the numerically observed distributions of
$\chi_{\ell}(\epsilon)$ for several $\ell$ at the center of the spectrum for
$U=1$ and $W=4.8$ (the critical point for the chain Hamiltonian (10a) - see
Appendix B), along with the results of carrying out this fit. In all cases,
the best-fit distribution is a uniform distribution on $[0,1]$. Though we only
show results for the center of the spectrum for $U=1$, we find similar
distributions at the transition throughout the spectrum, as well for
$U=2,3,4,5$. This observation that the amplitudes $\chi_{\ell}(\epsilon)$ are
uniformly distributed at the crossover will be key in the subsequent sections
when we derive expressions for the distribution of $g_{c}$ at the critical
point.

As a final note, the behavior of the distributions observed in Figures 6 and 7
is generic throughout the spectrum for $2\leq\ell\lessapprox 10$ (though the
peak of the distributions at 0 in the MBL phase is less pronounced for
$\ell\gtrapprox 6$). For the largest and smallest values of $\ell$, we do see
(not shown) deviation from the behavior in Fig. 6; we expect this to be caused
chiefly by finite-size and finite numerical precision effects. Going forward,
we assume that the qualitative behavior observed in Figures 6 and 7 describes
the asymptotic ($L\to\infty$) behavior of the system, and we attempt to
characterize our finite-size numerics by this expected behavior.

(a)

(b)

(c)

(d)

(e)

Figure 6:  Distributions of $\chi_{4}(\epsilon)$ (with respect to disorder
realizations) at $\epsilon=0.5$ for a variety of $W$, for a) $U=1$, b) $U=2$,
c) $U=3$, d) $U=4$, e) $U=5$. Each line is the observed probability
distribution $p(\chi_{4})$ (smoothed using a Gaussian kernel density
estimator) for a given disorder strength, with the color of the line
indicating the disorder strength. The dashed line traces out a uniform
distribution for comparison. We see a qualitative difference between the
distributions as we tune from the ergodic phase (dark blue) to deep in the MBL
phase (dark red); the distributions cross over from peaked at $1$ to peaked
near zero, with an intermediate critical regime of being approximately
uniform. All quantities are computed with $L=11$, $L_{b}=3$, in the $N=3$
sector, and the distributions are over 500 disorder realizations.

(a)

(b)

(c)

(d)

Figure 7: Histogram of amplitude ratios
$\chi_{\ell}(\epsilon)=\mathcal{A}_{\ell}(\epsilon)/\mathcal{A}_{\ell-1}(\epsilon)$
for energy density $\epsilon=0.5$ and interaction/disorder strength $U=1$,
$W=4.8$. Shown are the histograms for a) $\ell=4$, b) $\ell=5$, c) $\ell=6$,
d) $\ell=7$. Blue lines overlaid are the “best fit” distributions minimizing
the Akaike information (see text). The best fit is a uniform distribution in
all cases, suggesting $\chi_{i}\sim\text{Uni}[0,1]$ at the finite-size
crossover. All quantities are computed with $L=11$, $L_{b}=3$, in the $N=3$
sector, and the distributions are over 500 disorder realizations.

###  III.5 Connecting $g_{c}$ to $\mathcal{A}_{\ell}$

In section III.2, we found evidence that the non-hermitian many-body
delocalization transition is consistent with a non-hermitian avalanche
mechanism, and that $g_{c}\approx 1/\lambda-1/\lambda_{c}$ measured the
location of this transition, for the toy Hamiltonian (8). With the mapping of
the previous two sections, we are now in position to test if an analogous
relation holds for more realistic Hamiltonians (like the Hamiltonian (1)).

To that end, we take the Hamiltonian specified by (10b), add an imaginary flux
$i(L+1)g$, and compute $g_{c}(\epsilon)$ for a variety of energy densities
(see Appendix A for details on computing $g_{c}$). We do so for the same
disorder realizations as shown in the previous section, so that we have both
$g_{c}(\epsilon)$ and $\mathcal{A}_{\ell}(\epsilon)$ for the same disorder
realizations.

To connect $g_{c}(\epsilon)$ to the $\mathcal{A}_{\ell}(\epsilon)$, let us
consider an avalanche proliferating in this chain. As remarked in the last
section, the deterministic decay of hopping amplitudes $e^{-1/\lambda}$ in the
toy model (8) must be replaced by the random variables
$\mathcal{A}_{\ell}(\epsilon)$. Because of this, the ratio of matrix element
to bath level spacing (in analogy with (6)) is not necessarily monotonic, and
we must examine this ratio for every l-bit to determine delocalization. Due to
the gauge freedom, each of these bonds should also be enhanced by a factor of
$e^{(L+1)g}$ (since we can put all influence of $g$ onto any single bond in
the original Hamiltonian (10b), including the ones in $H_{\text{bc}}$). Thus,
the analogous condition to (6) (at energy density $\epsilon$) for all l-bits
should read:

$g\gtrapprox-\frac{1}{L+1}\min_{1\leq\ell\leq
L}\left[\log\mathcal{A}_{\ell}(\epsilon)+(\ell-1)\frac{\log 2}{2}\right]$ (14)

In words, the flux $i(L+1)g$ must make up for the maximal “difference” between
the (random) decay of the hoppings and the level spacing.

Based on the results of Section III.2, we expect that the exceptional points
$g_{c}(\epsilon)$ saturate this bound - that is, $g_{c}$ is exactly enough to
cover all of these “differences”. We test this criterion numerically by
computing

$\mathscr{A}(\epsilon):=-\min_{1\leq\ell\leq
L}\left[\log\mathcal{A}_{\ell}(\epsilon)+(\ell-1)\frac{\log 2}{2}\right]$ (15)

and $g_{c}(\epsilon)$ for each disorder realization, and comparing the two.
The result is very noisy, but there does appear to be a linear relationship
between these quantities - see Fig. 8 for some example scatter plots. We
observe a much clearer trend by averaging $g_{c}(\epsilon)$ and
$\mathscr{A}(\epsilon)$ over disorder realizations for multiple disorder
strengths $W$ \- see Fig. 9 for examples at the center of the spectrum. We are
able to fit a line to the scatter plot of $\left\langle
g_{c}(\epsilon)\right\rangle(W)$ vs
$\left\langle\mathscr{A}(\epsilon)\right\rangle(W)$ (where
$\left\langle\cdot\right\rangle$ denotes a disorder average) - these fits are
shown in Fig. 9 (and overlaid in Fig. 8), and Table 1 shows the resulting fit
parameters (with uncertainties 333Uncertainties on the fit parameters were
found via bootstrapping: for each disorder strength, we sample with
replacement from the disorder realizations, average $g_{c}(\epsilon)$ and
$\mathscr{A}(\epsilon)$, and perform a linear fit. The uncertainties on the
slope and intercept are the the standard deviation of each quantity over 1000
such resamplings.). We find that the disorder-averaged data is well-fit by a
line, and the slope agrees within error bars (at least at the center of the
spectrum) with the expected value of $1/12=.0825$. Curiously, the intercepts
of these fits are non-zero. At present, we do not have an explanation for
this, but we note in passing that the intercepts are roughly consistent with
$\log 2\approx 0.693$, which would correspond to a shift of about $\log 2$ per
non-bath site to $\mathscr{A}(\epsilon)$. In any case, there is unmistakably a
linear relationship between $\left\langle g_{c}\right\rangle$ and
$\left\langle\mathscr{A}(\epsilon)\right\rangle$, with (roughly) the correct
slope of $1/(L+1)$. This results support our conjecture in Section III.2, and
the idea that the location of $g_{c}$ is given by saturating the avalanche
criterion (14). This in turn means that the localization length $\xi=1/g_{c}$
captures the difference between the hopping decay and the level spacing.

In taking (15) as the “avalanche condition” determining $g_{c}$, we are
assuming that delocalization occurrs when all l-bits are delocalized by
hybridizing with the bath. It is a worthwhile question to ask whether this is
strictly necessary, and if delocalization of the system (indicated here by the
eigenvalues acuqiring a non-zero imaginary part) only requires a fraction of
the l-bits to hybridize with the bath. Requiring, say, only $L^{\prime}$ out
of $L$ l-bits to hybridize would entail replacing the minimum in (15) by the
function that selects the $(L^{\prime})$’th largest of the collection of $L$
values. In the extreme case of $L^{\prime}=1$, the minimum in (15) should be
replaced by a maximum; we find (not shown) that scatter plots analogous to
Fig. 8 are very noisy and uncorrelated, and disorder-averaged fits analogous
to those in Fig. 9 are very poor. Consequently, this extreme case is unlikely
to capture the physics (unlike the minimum, corresponding to $L^{\prime}=L$).

It is not clear whether there even should be a universal value of $L^{\prime}$
for all disorder realizations. If $L^{\prime}$ varies from realization to
realization, this could explain the non-zero intercepts in the linear fits of
Fig. 9; the intercept could capture the (disorder-averaged) variation between
the minimum and $(L^{\prime})$’th largest value of
$\mathcal{A}_{\ell}(\epsilon)+(\ell-1)\log 2/2$ across disorder realizations.
To attempt to address this question, we repeat the analysis presented above
with an average over l-bits, instead of a minimum - see Appendix G. There, we
find comparable results to those we have presented in this section, and the
linear fits generically have smaller intercepts. As averaging over l-bits does
not fix a particular universal value of $L^{\prime}$, this suggests that
$L^{\prime}$ may indeed not have a universal value across disorder
realizations. In spite of this - as we shall see in the next section - the
avalanche condition (15) (i.e., taking $L^{\prime}=L$) appears to be the most
successful in capturing the physics of $g_{c}$ in the original Hamiltonian (1)
without a bath, so we primarily focus on this case moving forward.

(a)

(b)

(c)

(d)

(e)

Figure 8: Example scatter plots of exceptional points $g_{c}(\epsilon)$ and
“avalanche parameter” $\mathscr{A}(\epsilon)$ (computed from the
$\mathcal{A}_{\ell}(\epsilon)$ found in Section III.4), slightly above the
critical point. Each point represents a single disorder realization (for a
fixed disorder strength), with the ordinate and abscissa $g_{c}$ and
$\mathscr{A}$ at the center of the spectrum ($\epsilon=0.5$), respectively.
The various panels show interaction/disorder strengths a) $U=1,W=5.5$, b)
$U=2,W=6.0$, c) $U=3,W=7.0$, d) $U=4,W=7.0$, e) $U=5,W=8.0$. We can see a
roughly linear relationship between $g_{c}(\epsilon)$ and
$\mathscr{A}(\epsilon)$, though there is significant noise. We quantify this
relationship by computing the statistical correlation
$\rho(X,Y)=\text{cov}(X,Y)/\sigma_{X}\sigma_{Y}$ (where cov is the covariance)
between $g_{c}(\epsilon)$ and $\mathscr{A}(\epsilon)$ \- this is marked in the
upper left of each plot. We see in all cases weak positive correlation
($\rho\sim 0.3$), consistent with the noisy linear relationship described
above. For reference, we also overlay the disorder averaged fits of Fig. 9
(dashed lines); these lines appear to roughly capture the relationship between
$g_{c}(\epsilon)$ and $\mathscr{A}(\epsilon)$. All data was taken from the
$N=3$ sector of a chain with $L=11$, $L_{b}=3$, and computed for 500 disorder
realizations.

(a)

(b)

(c)

(d)

(e)

Figure 9: Scatter plots of disorder-averaged exceptional points $g_{c}(\epsilon)$ and “avalanche parameter” $\mathscr{A}(\epsilon)$(computed from the $\mathcal{A}_{\ell}(\epsilon)$ found in Section III.4) for a variety of disorder strengths. Each point in the scatter plot corresponds to a specific disorder strength $W$ (the coloration indicating the numerical value of $W$), with the ordinate and abscissa the average (over disorder realizations) of $g_{c}$ and $\mathscr{A}$ at the center of the spectrum ($\epsilon=0.5$), respectively. The various panels show interaction strengths a) $U=1$, b) $U=2$, c) $U=3$, d) $U=4$, e) $U=5$. Dashed lines indicate a linear fit performed via orthogonal distance regression. We see good agreement between the fits and the data points, providing evidence that the real-complex/localization-delocalization transition is driven by non-hermitian avalanches. All data was taken from the $N=3$ sector of a chain with $L=11$, $L_{b}=3$. Averages were performed over 500 disorder realizations. | $U=1$ | $U=2$ | $U=3$ | $U=4$ | $U=5$
---|---|---|---|---|---
$\epsilon$ | Slope | Intercept | Slope | Intercept | Slope | Intercept | Slope | Intercept | Slope | Intercept
$0.2$ | $.098\pm.002$ | $-.895\pm.033$ | $.086\pm.002$ | $-.651\pm.031$ | $.091\pm.004$ | $-.737\pm.075$ | $.085\pm.003$ | $-.625\pm.044$ | $.093\pm.004$ | $-.694\pm.068$
$0.3$ | $.086\pm.001$ | $-.711\pm.023$ | $.080\pm.001$ | $-.586\pm.020$ | $.087\pm.003$ | $-.691\pm.052$ | $.083\pm.002$ | $-.619\pm.033$ | $.093\pm.003$ | $-.751\pm.048$
$0.4$ | $.084\pm.001$ | $-.683\pm.019$ | $.079\pm.001$ | $-.573\pm.018$ | $.086\pm.003$ | $-.664\pm.044$ | $.083\pm.002$ | $-.627\pm.032$ | $.094\pm.005$ | $-.769\pm.072$
$0.5$ | $.082\pm.001$ | $-.643\pm.018$ | $.079\pm.001$ | $-.574\pm.017$ | $.081\pm.002$ | $-.584\pm.038$ | $.083\pm.002$ | $-.613\pm.031$ | $.084\pm.004$ | $-.599\pm.066$
$0.6$ | $.082\pm.001$ | $-.643\pm.020$ | $.081\pm.001$ | $-.581\pm.018$ | $.089\pm.003$ | $-.703\pm.050$ | $.087\pm.002$ | $-.653\pm.037$ | $.100\pm.006$ | $-.798\pm.097$
$0.7$ | $.088\pm.002$ | $-.724\pm.025$ | $.082\pm.001$ | $-.586\pm.021$ | $.091\pm.004$ | $-.705\pm.059$ | $.091\pm.003$ | $-.663\pm.046$ | $.106\pm.004$ | $-.890\pm.070$
$0.8$ | $.092\pm.002$ | $-.775\pm.035$ | $.088\pm.002$ | $-.640\pm.030$ | $.091\pm.006$ | $-.651\pm.094$ | $.096\pm.004$ | $-.684\pm.073$ | $.103\pm.006$ | $-.748\pm.105$
Table 1: Best fit parameters for a linear fit of $\left\langle
g_{c}(\epsilon)\right\rangle$ vs
$\left\langle\mathscr{A}(\epsilon)\right\rangle$ (see the main text and Fig. 9
for details regarding the fits), across the spectrum for various interaction
strengths. The linear fits are performed using orthogonal distance regression,
using the standard errors of $\left\langle g_{c}(\epsilon)\right\rangle$ and
$\left\langle\mathscr{A}(\epsilon)\right\rangle$ as uncertainties. We estimate
the uncertainties on the fit parameters as the standard deviations on the fit
parameters across 1000 bootstrap resamplings. At weaker interaction strengths
and near the center of the spectrum, we see good agreement with the expected
slope of 1/12, while the parameters deviate from this expectation for stronger
interactions and towards the band edges. All data was taken from the $N=3$
sector of a chain with $L=11$, $L_{b}=3$. Averages were performed over 500
disorder realizations.

##  IV Distributions at the MBL crossover

We have sucessfully connected the localization length $\xi=1/g_{c}$ to the
length scale appearing in the avalanche model. We now seek to leverage this
relationship to predict the distribution of $g_{c}$ (with respect to disorder)
- the distribution of $\xi$ then follows as the inverse distribution, per the
relationship $\xi=1/g_{c}$. We focus here on the distribution at the MBL-
thermal crossover, as we have an analytic form for the distribution of the
amplitudes $\mathcal{A}_{\ell}(\epsilon)$ (based on the observations of
Section III.4).

Based on the results of Section III.5, we consider the following generalized
relation between $g_{c}$ and the $\mathcal{A}_{\ell}$’s:

$g_{c}(\epsilon)\approx\frac{1}{F}\max_{1\leq\ell\leq
N_{\ell}}\left[-\log\mathcal{A}_{\ell}(\epsilon)-(\ell-1)\alpha\right]+g_{0}$
(16)

where $iFg$ is the total flux through the system, $N_{\ell}$ the number of
l-bits in the system, and $\alpha=\log d/2$, where $d$ is the on-site Hilbert
space dimension (in fermionic models, $d=2$). We also allow for an offset
$g_{0}$, based on the results of Section III.5. Note that we have written the
relation in terms of a maximum here, to ease analytic computation of the
distribution.

Let us assume this relation for $g_{c}$ holds for each disorder realization
(i.e., we neglect the noise observed in Fig. 8), and let
$P_{F}^{N_{\ell}}(g_{c}\leq x\,|\,g_{0})$ denote the cumulative density
function for $g_{c}$ given by this relation. As observed in Section III.4, we
have $\mathcal{A}_{\ell}(\epsilon)=\prod_{i=1}^{\ell}\chi_{i}(\epsilon)$ with
$\chi_{i}(\epsilon)$ uniformly distributed on $[0,1]$ at the MBL crossover.
Using this, we arrive at the following form for the cdf

$\displaystyle P_{F}^{N_{\ell}}(g_{c}\leq
x\,|\,g_{0})=1-e^{-\tilde{x}}\left[1+\tilde{x}\sum_{j=1}^{N_{\ell}-1}\frac{e^{-j\alpha}}{j!}(\tilde{x}+j\alpha)^{j}\right]$

(17)

where $\tilde{x}=F(x-g_{0})$ \- see Appendix E for details on the derivation
of this result. Crucially, we note here that the derivation in Appendix E
assumes the $\chi_{i}$ are independently and identically distributed. This is
not strictly true, as the distributions of $\chi_{i}$ do vary in $i$ (as
remarked in Section III.4), and all the $\chi_{i}(\epsilon)$ arise from the
same disorder realization, and are thus almost certainly correlated. We
neglect these two facts here.

We now would like to see if the analytic distribution (17) describes the
distribution of $g_{c}$’s (and thereby $\xi$’s) at the MBL crossover in the
original model, (1). To that end, we compute the histogram of
$g_{c}(\epsilon)$ at the (energy-resolved) critical disorder strength
identified from finite-size scaling (see Appendix B). We cannot compare the
numerically observed distributions directly to (17) however, as the
Hamiltonian (1) doesn’t explicitly include a thermal bath/grain. The size of
the bath is thus itself also a random variable; we account for this by
modifying (17) to a mixture distribution of the form

$P_{\bm{\beta},g_{0}}(g_{c}\leq
x)=\sum_{i=1}^{L-1}\beta_{i}P_{F}^{L-i}(g_{c}\leq x\,|\,g_{0})$ (18)

with $0\leq\beta_{i}\leq 1$, $\sum_{i}\beta_{i}=1$. In this distribution,
$\beta_{i}$ is the probability the bath/thermal grain is of size $i$, and
$P_{F}^{L-i}$ is the CDF of $g_{c}$ conditioned on the bath being size $i$.
Such a distribution follows from the law of total probability applied to the
CDF of $g_{c}$.

We fit a mixture distribution of the form (18) to the numerically observed
CDF. We do so by computing a kernel density estimator $\hat{f}$ of the
numerically observed CDF (using a Gaussian kernel), sampling this estimator at
various points $\\{x_{i}\\}_{i=1}^{n}$, and numerically minimizing the squared
error

$\delta(\bm{\beta},g_{0})=\sum_{i=1}^{n}\left(P_{\bm{\beta},g_{0}}(g_{c}\leq
x_{i})-\hat{f}(x_{i})\right)^{2}$ (19)

subject to the constraints $0\leq\beta_{i}\leq 1$, $\sum_{i}\beta_{i}=1$. The
results of this fit at the center of the spectrum are shown in Fig 10. We see
good overall agreement between the observed CDF and the fitted mixture
distribution, though the fitted distributions have longer tails than the
observed distributions (as evidenced by the slower approach of the CDF to 1 -
see also Appendix F). This suggests that the distributions of $g_{c}$ at the
critical point are well-described by an avalanche-like picture of
delocalization (which gives rise to the distribution (18)). The distribution
of $\xi$, obtained by computing the inverse distribution of (18), thus
contains information about the distribution of the hopping amplitudes
$\mathcal{A}_{\ell}$.

Curiously, we find all the fits have $g_{0}\approx 0$; the offset does not
appear to play a role in the physics without the bath. In light of the
discussion at the end of the previous section, the lack of offset would
suggest that all l-bits ($L^{\prime}=L$) must hybridize when there is no
explicit bath in the system. In Appendix G, we derive an analytic distribution
of $g_{c}$ assuming an average avalanche condition and perform similar fits;
the agreement between fitted and numerically observed distributions is not as
good for the average avalanche condition (especially at larger interaction
strengths). This further suggests that the physics in the system without a
bath is captured by requiring all l-bits to hybridize ($L^{\prime}=L$).

As a final remark, we note that the fit quality of the mixture distribution
(18) is even better towards the band edges - see Appendix F for examples and
discussion.

(a)

(b)

(c)

(d)

(e)

Figure 10: Numerically observed cumulative density functions of $g_{c}$ for
$\epsilon=0.5$ at the critical point, for interaction strengths a) $U=1$
($W=4.8$) , b) $U=2$ ($W=5.4$), c) $U=3$ ($W=5.9$), d) $U=4$ ($W=6.1$), e)
$U=5$ ($W=6.7$). The blue solid lines overlaid are best fit mixture
distributions of the form (18) (see main text for details). We see good
agreement between fit and numerical observation, suggesting that our
distribution derived from a non-hermitian avalanche criterion accurately
describes the transition point. The parameters extracted from the fit are a)
$\beta_{5}=0.366$, $\beta_{6}=0.634$ b) $\beta_{5}=0.769$, $\beta_{6}=0.231$,
c) $\beta_{6}=0.435$, $\beta_{7}=0.565$, d) $\beta_{4}=0.827$,
$\beta_{5}=0.173$ e) $\beta_{3}=0.243$, $\beta_{4}=0.757$, and $g_{0}=0$ in
all cases. All quantities are computed in chains of size $L=12$ in the half-
filling sector, and from 10000 disorder realizations. When computing the
squared error (19), we sample $n=100$ points uniformly spaced between the
smallest and largest observed values of $g_{c}$.

##  V Conclusion & Discussion

We have shown how to connect a localization length $\xi$, defined by
introducing an imaginary vector potential to chain of interacting disordered
fermions, to the length scale $\lambda$ characterizing the decay of matrix
elements in the avalanche model of delocalization. We have also derived an
analytic form of how the $g_{c}$’s are distributed at the MBL transition in
finite systems, and shown that this form describes well the observed
histograms. The distributions of $\xi$ at the transition, which can be derived
from those of $g_{c}$, thus contain information about the distribution of
avalanche-like hoppings, giving insight into the avalanche mechanism at the
transition.

Determining $\xi$, while computationally difficult, has the advantage of not
requiring one to explicitly construct the l-bits. That is, a given $H$ with a
fixed disorder realization has a well-defined $g_{c}$, and thereby $\xi$, for
each energy density $\epsilon$. Though our numerical verification of $g_{c}$
saturating the avalanche condition (14) made use of a particular construction
of the l-bits, the agreement of the numerically observed histograms to our
fitted CDFs (18) suggests the form of avalanche criteria for $g_{c}$ does not
depend on how the l-bits are constructed. This method thus allows one to
obtain information about the localization of the underlying $g=0$ Hamiltonian
whilst sidestepping issues with defining the l-bits.

The connection of $\xi$ to the localization properties of the underlying $g=0$
system relies on the proliferation of avalanches enhanced by the imaginary
flux through the system. It is therefore an intereresting question to ask how
a $\xi$ defined analogously for “clean” systems differs, as there are no rare
regions driving thermalization in these systems. Recent works Villalonga and
Clark (2020a, b); Abanin et al. (2021); Crowley and Chandran (2022);
Morningstar et al. (2022) have also pointed out the difficulties of analyzing
MBL in numerically accessible systems, and how the finite-size MBL-thermal
crossover is driven by the proliferation of so-called “many-body resonances”,
in which two (or more) eigenstates hybridize under a local perturbation,
rather than avalanches. A further analysis of how $g_{c}$ is affected by the
presence of these resonances - e.g. how it behaves in the resonance model of
Ref Crowley and Chandran (2022) \- could provide additional insight, and
quantify to what extent the length scale $\xi$ suffers from the finite-size
effects plaguing many MBL systems.

###### Acknowledgements.

LO thanks Christopher David White for his guidance and advice during the early
stages of this project, as well as Dan Borgnia for useful discussions. GR is
grateful for support from the Simons Foundation as well as support from the
NSF DMR grant number 1839271, and from the IQIM, an NSF Physics Frontiers
Center. This work was performed in part at Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-1607611.

##  Appendix A Numerical Procedure for Obtaining $g_{c}$

Here, we briefly summarize our numerical procedure for obtaining the
exceptional points $g_{c}$.

For a given instance of a Hamiltonian (i.e., a fixed disorder realization),
the exceptional point $g_{c}(\epsilon)$ at energy density $\epsilon$ is
defined as:

$g_{c}(\epsilon)=\inf\left\\{g\,|\,g>0,\,\text{Im}[E({\epsilon})]\neq
0\right\\}$ (20)

where $E({\epsilon})$ is the eigenvalue at energy density $\epsilon$. The
length scale $\xi(\epsilon)$ is then defined by
$\xi(\epsilon)=1/g_{c}(\epsilon)$ \- from here on we refer to $\xi(\epsilon)$
as the localization length.

For finite-size systems, the energies are discrete, and we must consider the
eigenvalue closest to energy density $\epsilon$, given by

$E({\epsilon})=\arg\min_{E\in\mathcal{S}\\{H\\}}\left\lvert\frac{\text{Re}\left[E\right]-\text{Re}\left[E_{\min}\right]}{\text{Re}\left[E_{\max}\right]-\text{Re}\left[E_{\min}\right]}-\epsilon\right\rvert$
(21)

where $\mathcal{S}\\{H\\}$ is the spectrum of $H$, and $E_{\max}$, $E_{\min}$
are the eigenvalues with the maximal and minimal real parts (respectively).

To find $g_{c}(\epsilon)$ numerically for a fixed disorder realization, we
repeatedly increment $g$ by a fixed step $\delta g$ (which we take to be .01),
exactly diagonalizing the Hamiltonian at each value of $g$, until
$|\text{Im}\left[E(\epsilon)\right]|$ is larger than some tolerance (which we
take to be $10^{-8}$). Note this numerical approach differs from previous
studies that examine all eigenvalues within some energy window to determine
$g_{c}$. We choose to sample only one state per disorder realization to
sidestep subtleties involving correlations of eigenvectors from the same
disorder distributions Rodriguez et al. (2011). We have verified that this
approach gives the same disorder-averaged results as averaging over an energy
window.

##  Appendix B Identifying the MBL transition via finite-size scaling

In this Appendix, we show finite-size scaling collapse on the disorder-
averaged localization lengths to identify the MBL transition for the
Hamiltonian (1). We find results in good agreement with those of Heußen et al.
(2021).

To that end, we compute $g_{c}(\epsilon)$ (as described in Appendix A) for the
Hamiltonian (1), with bare hopping $t=1$ (making it the energy scale of our
system) and on-site disorders sampled from a uniform distribution
$w_{i}\sim\text{Uni}[-W,W]$. Figure 11 shows finite size-scaling of
$\bar{\xi}/L=\left\langle g_{c}\right\rangle^{-1}/L$ (where
$\left\langle\cdot\right\rangle$ denotes a disorder average) at the center of
the spectrum ($\epsilon=0.5$) in the half-filling sector. We examine a variety
of interaction strengths $U=1,\,2,\,3,\,4,\,5$, at system sizes
$L=8,\,10,\,12$, and average over 200 disorder realizations. We obtain good
collapse of the data to the form

$\frac{\bar{\xi}}{L}=f\left((W-W_{c})L^{1/\nu}\right)$ (22)

and in all cases find critical exponents $\nu\approx 1$. The resulting phase
portraits in the $\epsilon$-$W$ plane are shown in Fig. 12, with the
transition points identified from the finite size scaling marked. We see the
transition points capture the change from large (order of system size) to
small localization length as we tune the disorder strength. The shape of the
mobillity edge in the phase diagram also mirrors that of other finite size
studies Luitz et al. (2015), and the critical disorder strengths agree well
with Heußen et al. (2021) .

(a) $U=1$

(b) $U=2$

(c) $U=3$

(d) $U=4$

(e) $U=5$

Figure 11: Finite-size scaling collapse for $\bar{\xi}/L=\left\langle
g_{c}\right\rangle^{-1}/L$ at energy density $\epsilon=0.5$ for a) $U=1$ , b)
$U=2$, c)$U=3$, d) $U=4$, e) $U=5$. The critical exponent extracted in all
cases is $\nu\approx 1$; the critical disorder strengths are $W_{c}\approx
4.8,5.4,5.9,6.1,6.7$ for $U=1,2,3,4,5$, respectively. All quantities are
computed in the half-filling sector, and averaged over 200 disorder
realizations.

(a) $U=1$

(b) $U=2$

(c) $U=3$

(d) $U=4$

(e) $U=5$

Figure 12: Phase portraits for $\bar{\xi}/L=\left\langle
g_{c}\right\rangle^{-1}/L$ in the $\epsilon-W$ plane, for a) $U=1$ , b) $U=2$,
c) $U=3$, d) $U=4$, e) $U=5$, for a chain of size $L=12$. White triangles mark
the critical disorder strengths obtained from our finite size scaling
collapse. All quantities are computed in the half-filling sector, and averaged
over 200 disorder realizations.

##  Appendix C Additional Results Regarding Non-Hermitian Avalanches

Here, we show additional numerical results supporting the existence of the
non-hermitian avalanche mechanism, as in Fig. 3. Results for a matrix element
decay of $e^{-1/\lambda}=0.2$ and $e^{-1/\lambda}=0.1$ are shown in Fig. 13
and Fig. 14, respectively. We see crossovers analogous to those observed in
the main text, again occurring at the expected value of
$g=1/\lambda-1/\lambda_{c}$.

(a)

(b)

(c)

Figure 13: Disorder-averaged localization metrics of the eigenstates of the
non-hermitian avalanche Hamiltonian (8) for $e^{-1/\lambda}=0.2$ and various
$L$, as a function of $g$. Going clockwise: a) half-chain entanglement entropy
$\langle S_{E}\rangle$, b) Inverse participation ratio $\langle IPR\rangle$,
and c) Fraction $f_{\text{comp}}$ of disorder realizations whose central
eigenvalue ($\epsilon=0.5$) is non-real. As in Fig. 3, we see a crossover from
localized to delocalized behavior in panels a) and b) at the expected value of
$g$ (dashed lines), along with a corresponding real-complex crossover in panel
c).
All quantities were computed exactly as in Fig. 3.

(a)

(b)

(c)

Figure 14: Disorder-averaged localization metrics of the eigenstates of the
non-hermitian avalanche Hamiltonian (8), for $e^{-1/\lambda}=0.1$ and various
$L$, as a function of $g$. Going clockwise: a) half-chain entanglement entropy
$\langle S_{E}\rangle$, b) Inverse participation ratio $\langle IPR\rangle$,
and c) Fraction $f_{\text{comp}}$ of disorder realizations whose central
eigenvalue ($\epsilon=0.5$) is non-real for various $L$ and $g$. We observe
similar crossovers as seen in Fig. 3 and Fig. 13.

##  Appendix D Diagonalizing many-body hamiltonians via generalized
displacement transformations

In this appendix, we detail a diagonalization algorithm that iteratively
applies a series of transformations eliminating “off-diagonal” terms from a
given fermionic Hamiltonian.

Consider fermions hopping on a 1-dimensional lattice with sites
$S=\\{i\\}_{i=1}^{L}$. Let $c_{i}$, $c_{i}^{{\dagger}}$, and $n_{i}$ be the
fermionic annhiliation, creation, and number operators at a site $i$,
respectively. A general (hermitian) Hamiltonian then has the form:

$H=V+\sum_{A,B}t_{AB}I_{AB}$ (23)

where $V$ is a diagonal operator, $A,B$ are disjoint subsets of $S$, $t_{AB}$
is a diagonal operator acting only on sites in $S\setminus A\cup B$, and
$I_{AB}$ is the operator for hopping between $A$ and $B$:

$I_{AB}=\left(\prod_{a\in A}c_{a}^{{\dagger}}\right)\left(\prod_{b\in
B}c_{b}\right)+\text{ h.c. }$ (24)

We neglect here the possibility of complex phases multiplying the hopping
terms - our algorithm presented below is straightforwardly modified in this
case, and the essential results do not change.

Suppose we now want to eliminate the hopping $I_{AB}$ via a unitary change of
basis. We can do this by the following transformation

$H\mapsto e^{\mathcal{O}_{AB}J_{AB}}\,H\,e^{-\mathcal{O}_{AB}J_{AB}}$ (25)

where

$J_{AB}=\left(\prod_{a\in A}c_{a}^{{\dagger}}\right)\left(\prod_{b\in
B}c_{b}\right)-\text{h.c.}$ (26)

and $\mathcal{O}_{AB}$ is a diagonal operator acting only on sites in
$S\setminus A\cup B$ satisfying

$\tan\left(2\mathcal{O}_{AB}\right)=\frac{2t_{AB}}{\Delta_{AB}V}$ (27)

where

$\Delta_{AB}V:=V\big{\rvert}_{\begin{subarray}{c}A\text{ occupied},\\\ B\text{
unoccupied}\end{subarray}}-V\big{\rvert}_{\begin{subarray}{c}B\text{
occupied},\\\ A\text{ unoccupied}\end{subarray}}$ (28)

Application of the transformation (25) to the Hamiltonian (23) results in a
Hamiltonian with no terms of the form $\mathcal{D}\,I_{AB}$ present (where
$\mathcal{D}$ is a diagonal operator). Generically, such a Hamiltonian will
include new diagonal terms, as well as new hopping terms (created by the
action of the transformation on the other hopping terms in (23)).

Now, we can repeatedly apply transformations of the form (25) to eliminate all
hopping terms in the original Hamiltonian (23). As mentioned previously, each
application of such a transformation generates new hopping terms; in the case
of repeated transformations, we can actually re-introduce hoppings that were
eliminated by previous transformations. This is not a problem, as the terms
are always re-introduced with a smaller “magnitude”. To be more precise, let
us write out the unitary map in (25) in closed form:

$e^{\pm\mathcal{O}_{AB}J_{AB}}=\mathbb{I}+[\cos(\mathcal{O}_{AB})-\mathbb{I}]P_{AB}\pm\sin(\mathcal{O}_{AB})J_{AB}$

(29)

where $P_{AB}$ is the projector onto the subspace in which all sites in either
$A$ or $B$ (but not both) are filled:

$\displaystyle P_{AB}=\left(\prod_{a\in A}n_{a}\right)\left(\prod_{b\in
B}(\mathbb{I}-n_{b})\right)+\left(\prod_{b\in B}n_{b}\right)\left(\prod_{a\in
A}(\mathbb{I}-n_{a})\right)$

(30)

In the form (29), we can see that every new hopping term introduced by a
transformation of the form (25) will include an operator multiplication by a
trigonemtric function of $\mathcal{O}_{AB}$, which will have operator norm
$<1$. If we iterate such transformations, we can see that the re-introduced
hopping terms will converge (in operator norm) to zero. Thus, in the limit of
infinitely many transformations of the form (25) applied, all the hopping
(non-diagonal) terms will aproach zero, and we will be left with a diagonal
Hamiltonian. Numerically, we must cut the procedure off when the largest
(given by magnitude of largest matrix element) off-diagonal term has norm
below some threshold - we choose this to be $10^{-8}$.

Note also that by applying the same sequence of transformations to the
creation/annhilation operators, we can obtain creation/annhilation operators
for the local integrals of motion for the Hamiltonian (23) (assuming it is in
the MBL phase). Additionally, if we track the evolution of the hopping terms
as we iterate these transformations - or leave certain hopping terms untouched
by our transformations - we will see the amplitudes evolve as a function of
iteration number. This is reminiscent of a discretized version of the Wegner
flow employed by Quito et al. (2016). Indeed, this intuition is how we obtain
the “avalanche-like” hopping amplitudes described in the main text.

As a final note, this algorithm can be generalized to non-hermitian
Hamiltonians. In such a case, the transformation to eliminate a term of the
form

$\displaystyle I^{\alpha}_{AB}=\left(\prod_{a\in
A}c_{a}^{{\dagger}}\right)\left(\prod_{b\in
B}c_{b}\right)+\alpha\left(\prod_{b\in
B}c_{b}^{{\dagger}}\right)\left(\prod_{a\in A}c_{a}\right)$

(31)

\- where $\alpha$ is diagonal on $S\setminus A\cup B$ \- is constructed
analogously, but now with
$\cos(\mathcal{O}_{AB})\mapsto\cos(\sqrt{\alpha}\,\mathcal{O}_{AB})$ and
$\sin(\mathcal{O}_{AB})\mapsto\sin(\sqrt{\alpha}\,\mathcal{O}_{AB})/\sqrt{\alpha}$
in (27) and (29). Note that, if any of the eigenvalues of $\alpha$ are
negative, $\cos(\sqrt{\alpha}\,\mathcal{O}_{AB})$ and
$\sin(\sqrt{\alpha}\,\mathcal{O}_{AB})$ are no longer necessarily bounded in
operator norm, and convergence of the algorithm is no longer assured. This is
why it is crucial to use the gauge freedom to eliminate the non-hermiticity
from the part of the Hamiltonian we are diagonalizing, as described in the
main text.

These transformations (25), and the iterative procedure described here, are
examples of the “displacement transformations” and diagonalization algorithm
first proposed by Rademaker and Ortuño Rademaker and Ortuño (2016).

##  Appendix E Derivation of Cumulative Density Function at the MBL transition

In this appendix, we seek to derive an analytic expression for the
distribution of $g_{c}$ at the MBL transition, based on the (generalized)
avalanche condition (16). As it turns out, it is easier to derive the
cumulative density function (CDF) of $g_{c}$. To that end, let us consider the
following more general problem.

Let $\\{x_{i}\\}_{i=1}^{N}$ be a collection of random variables, and $\alpha$
a real number. Consider the following extremization problem:

$M_{L}=\max_{1\leq j\leq
L}\left\\{\sum_{i=1}^{j}x_{i}-(j-1)\alpha\right\\}:=\max_{1\leq j\leq L}m_{i}$
(32)

We wish to find the CDF of $M_{L}$, given the probability densities
$\rho_{i}(x_{i})$ of the $x_{i}$’s.

To write down the generic form of this CDF, we first make the following
observation:

$P(M_{L}\leq x)=P(m_{1}\leq x,m_{2}\leq x,\ldots,m_{L}\leq x)$ (33)

That is, $M_{L}\leq x$ iff $m_{j}\leq x$ for every value of $j$. Futhermore,
we can see that $m_{j}\leq x$ iff $x_{j}\leq
x+(j-1)\alpha-\sum_{i=1}^{j-1}x_{i}$. Hence:

$P(m_{j}\leq x)=P\left(x_{j}\leq x+(j-1)\alpha-\sum_{i=1}^{j-1}x_{i}\right)$
(34)

The right hand side is the CDF of $x_{j}$, evaluated at a point that depends
on the other $x_{i}$. Using this in the right side of (33), we obtain an
expression for $P(M_{L}\leq x)$ in terms of the PDFs of the $x_{j}$’s (i.e.,
the $\rho_{j}$’s):

$P(M_{L}\leq
x)=\prod_{j=1}^{L}\int_{-\infty}^{x+(j-1)\alpha-\sum_{i=1}^{j-1}x_{i}}\thinspace\mathrm{d}x_{j}\,\rho_{j}(x_{j})$
(35)

At this point, we can compute the integral over $x_{L}$ to obtain an integral
expression involving the CDF of $x_{L}$, but we cannot go any farther without
knowledge of the $\rho_{i}$’s.

To that end, let us now take each of the $x_{i}$’s to be exponentially
distributed - that is, $\rho_{i}(x_{i})=e^{-x_{i}}$ for $x\geq 0$, and $0$ for
$x<0$. We can perform these integrals for any finite value of $L$, and obtain:

$\displaystyle P(M_{L}\leq
x)=1-e^{-x}\left[1+x\sum_{j=1}^{L-1}\frac{e^{-j\alpha}}{j!}(x+j\alpha)^{j-1}\right]$

(36)

This is the desired CDF of $M_{L}$, assuming $x_{i}\sim\text{Exp}[1]$ for each
$i$. The CDF for a more general exponential distribution
$x_{i}\sim\text{Exp}[\lambda]$ can be obtained by making the replacements
$x,\,\mapsto\lambda x$, $\alpha\mapsto\lambda\alpha$.

Let us now connect back to the exceptional points. We have the generalized
“avalanche criterion” for $g_{c}$:

$g_{c}\approx\frac{1}{F}\max_{1\leq\ell\leq
N_{\ell}}\left[-\log\mathcal{A}_{\ell}(\epsilon)-(\ell-1)\alpha\right]+g_{0}$
(37)

Comparing (32) with this generalized criterion, and assuming
$\mathcal{A}_{\ell}(\epsilon)=\prod_{i=1}^{\ell}\chi_{i}(\epsilon)$, we see
that $g_{c}=\frac{1}{F}M_{N_{\ell}}$, with $x_{i}=-\log\chi_{i}(\epsilon)$.
Assuming $\chi_{i}\sim\text{Uni}[0,1]$ at the transition (as observed in
Section III.4), it follows that $-\log\chi_{i}(\epsilon)\sim\text{Exp}[1]$.
Therefore, the result (36) with the replacements $x\mapsto F(x-g_{0})$,
$L\mapsto N_{\ell}$ and $\alpha=\log 2/2$ should describe the distribution of
$g_{c}(\epsilon)$ at the MBL transition.

Note that the general formula (35), with appropriate replacements, gives the
distribution for $g_{c}$ whenever an avalanche criterion of the form (37)
holds. In particular, if the avalanche criteria holds in the MBL phase, this
formula should still give the distribution of $g_{c}$. Of course, the
$\chi_{i}(\epsilon)$ are not necessarily uniformly distributed away from the
transition; nonetheless, replacing $\rho_{i}$ by the appropriate distribution
of $-\log\chi_{i}(\epsilon)$ should yield the correct distribution for
$g_{c}(\epsilon)$.

##  Appendix F Additional results from Numerical Fits of Probability
distributions

In this appendix, we show some additional results from the numerical fits of
the CDFs of $g_{c}(\epsilon)$ described in section IV.

Figure 15 shows the numerically observed probability density function for
$g_{c}(\epsilon)$, along with the analytic PDF obtained by differentiating the
best fit mixture distribution of section IV, for the data in Fig 10. We can
see that the PDFs generated from our fit have thicker tails (i.e., approach
zero more slowly) than the numerically observed PDFs, and that the peak of the
numerically observed PDF is larger than the peak of the fitted distribution.
Despite these discrepancies, the PDFs from the fits still capture the observed
behavior qualitatively well, especially at small $g_{c}$ and for $U=3$.

(a)

(b)

(c)

(d)

(e)

Figure 15: Numerically observed probability density functions of $g_{c}$ for
$\epsilon=0.5$ at the critical point, for interaction strengths a) $U=1$
($W=4.8$) , b) $U=2$ ($W=5.4$), c) $U=3$ ($W=5.9$), d) $U=4$ ($W=6.1$), e)
$U=5$ ($W=6.7$). The blue solid lines overlaid are the pdfs obtained from
differentiating the best fit mixture distributions in Fig. 10.

As mentioned in the main text, we observed the numerical fits to be better
towards the edges of the spectrum. To exemplify this, Fig. 16 shows the fits
for the energy densities whose fitted mixture distribution had the lowest
squared error (defined by (19)), and Fig. 17 shows the corresponding PDFs. We
see excellent agreement between the observed and fitted distributions,
especially for $U\geq 2$.

(a)

(b)

(c)

(d)

(e)

Figure 16: Numerically observed CDFs at the critical point, with the best fit
mixture distribution overlaid, for a) $U=1$, $\epsilon=0.3$ ($W=3.4$), b)
$U=2$, $\epsilon=0.25$ ($W=3.9$), c) $U=3$, $\epsilon=0.8$ ($W=3.0$), d)
$U=4$, $\epsilon=0.3$ ($W=4.9$), e) $U=5$, $\epsilon=0.3$ ($W=5.7$). The
energy densities shown are the energy whose fitted distribution minimized the
squared error (19), for each interaction strength. We can see excellent
agreement between observation and fit, suggesting that the non-hermitian
avalanche mechanism developed in the main text is an even better description
of delocalization at the band edges. The parameters extracted from the fit are
a) $\beta_{8}=0.756$, $\beta_{9}=0.244$ b) $\beta_{7}=0.794$,
$\beta_{8}=0.206$, c) $\beta_{4}=0.481$, $\beta_{5}=0.518$, d)
$\beta_{5}=0.313$, $\beta_{6}=0.527$, $\beta_{8}=0.125$, $\beta_{9}=.035$, e)
$\beta_{1}=0.036$, $\beta_{4}=0.124$, $\beta_{5}=0.840$, and $g_{0}=0$ in all
cases. All quantities are computed exactly as in Fig. 10.

(a)

(b)

(c)

(d)

(e)

Figure 17: Numerically observed PDFs at the MBL transition, with the PDFs from
the best fit mixture distribution overlaid for, a) $U=1$, $\epsilon=0.3$
($W=3.4$), b) $U=2$, $\epsilon=0.25$ ($W=3.9$), c) $U=3$, $\epsilon=0.8$
($W=3.0$), d) $U=4$, $\epsilon=0.3$ ($W=4.9$), e) $U=5$, $\epsilon=0.3$
($W=5.7$). We see excellent agreement again, confirming what we observed in
Fig 16. In contrast with Fig. 15, we also see that the tails of the observed
distribution are well-described by the PDFs from the fits. All quantities are
computed exactly as in Fig. 10 and 16.

##  Appendix G Average Avalanche Condition

In the process of relating $g_{c}(\epsilon)$ to the amplitudes
$\mathcal{A}_{\ell}(\epsilon)$, we found empirically that $\left\langle
g_{c}(\epsilon)\right\rangle$ was also well-described by a linear relationship
with:

$\bar{\mathscr{A}}(\epsilon):=-\frac{1}{L}\sum_{\ell=1}^{L}\left[\log\mathcal{A}_{\ell}(\epsilon)+(\ell-1)\frac{\log
2}{2}\right]$ (38)

This quantity is an average over l-bits of the difference between hopping and
level spacing, instead of the extremum present in (15). As discussed in
Section III.5, such a relationship may be better suited to describing a
delocalization mechanism in which the number of l-bits that need to hybridize
with the bath for the system to undergo a real-complex transition varies
between disorder realizations. Figure 18 shows a scatterplot of
$g_{c}(\epsilon)$ versus $\bar{\mathscr{A}}(\epsilon)$ for individual disorder
realizations, akin to Fig. 8. As in Fig. 8, we see a noisy, approximately
linear relationship. We perform linear fits of $\left\langle
g_{c}(\epsilon)\right\rangle(W)$ vs
$\left\langle\bar{\mathscr{A}}(\epsilon)\right\rangle(W)$, in analogy with
Fig. 9; the results of these fits are shown in Fig. 19, and the extracted fit
parameters are shown in Table 2. The quality of fit is comparable to those
shown in the main text for $\left\langle g_{c}(\epsilon)\right\rangle$ vs
$\left\langle\mathscr{A}(\epsilon)\right\rangle$, though the fit parameters
are different.

(a)

(b)

(c)

(d)

(e)

Figure 18: Example scatter plots of exceptional points $g_{c}(\epsilon)$ and
“avalanche parameter” $\bar{\mathscr{A}}(\epsilon)$ (given by (38)), as in
Fig. 8. Each point represents a single disorder realization (for a fixed
disorder strength), with the ordinate and abscissa $g_{c}$ and
$\bar{\mathscr{A}}$ at the center of the spectrum ($\epsilon=0.5$),
respectively. The various panels show interaction/disorder strengths a)
$U=1,W=5.5$, b) $U=2,W=6.0$, c) $U=3,W=7.0$, d) $U=4,W=7.0$, e) $U=5,W=8.0$.
Similar to Fig. 8, we see a noisy, seemingly linear relationship between
$g_{c}(\epsilon)$ and $\mathscr{A}(\epsilon)$, with weak positive correlation
$\rho\sim 0.3$. For reference, we also overlay the disorder averaged fits of
Fig. 19 (dashed lines); these lines appear to roughly capture the relationship
between $g_{c}(\epsilon)$ and $\bar{\mathscr{A}}(\epsilon)$. All data was
computed exactly as in Fig. 8.

(a)

(b)

(c)

(d)

(e)

Figure 19: Scatter plot of disorder-averaged $g_{c}(\epsilon)$ and $\bar{\mathscr{A}}(\epsilon)$ for a variety of disorder strengths. The scatter plot is made exactly as in in Fig. 9 with $\bar{\mathscr{A}}$ replacing $\mathscr{A}$. Shown are results at the center of the spectrum ($\epsilon=0.5$) for interaction strengths a) $U=1$, b) $U=2$, c) $U=3$, d) $U=4$, e) $U=5$. Comparing with the fits in Fig. 9, we see comparable fit quality but different fit parameters. | $U=1$ | $U=2$ | $U=3$ | $U=4$ | $U=5$
---|---|---|---|---|---
$\epsilon$ | Slope | Intercept | Slope | Intercept | Slope | Intercept | Slope | Intercept | Slope | Intercept
$0.2$ | $.099\pm.002$ | $-.332\pm.018$ | $.093\pm.002$ | $-.224\pm.016$ | $.096\pm.004$ | $-.232\pm.041$ | $.094\pm.002$ | $-.198\pm.024$ | $.099\pm.004$ | $-.199\pm.041$
$0.3$ | $.094\pm.001$ | $-.315\pm.014$ | $.092\pm.001$ | $-.259\pm.012$ | $.101\pm.003$ | $-.326\pm.033$ | $.095\pm.002$ | $-.267\pm.020$ | $.103\pm.003$ | $-.318\pm.030$
$0.4$ | $.095\pm.001$ | $-.340\pm.012$ | $.095\pm.001$ | $-.296\pm.011$ | $.101\pm.003$ | $-.335\pm.026$ | $.100\pm.002$ | $-.315\pm.019$ | $.107\pm.004$ | $-.361\pm.041$
$0.5$ | $.093\pm.001$ | $-.324\pm.011$ | $.095\pm.001$ | $-.294\pm.011$ | $.100\pm.003$ | $-.323\pm.026$ | $.097\pm.002$ | $-.289\pm.019$ | $.096\pm.004$ | $-.246\pm.042$
$0.6$ | $.092\pm.001$ | $-.301\pm.012$ | $.097\pm.001$ | $-.298\pm.012$ | $.100\pm.003$ | $-.311\pm.029$ | $.097\pm.002$ | $-.263\pm.020$ | $.111\pm.006$ | $-.339\pm.058$
$0.7$ | $.097\pm.001$ | $-.323\pm.014$ | $.095\pm.001$ | $-.252\pm.014$ | $.098\pm.003$ | $-.254\pm.035$ | $.100\pm.003$ | $-.231\pm.029$ | $.111\pm.004$ | $-.323\pm.042$
$0.8$ | $.094\pm.002$ | $-.258\pm.019$ | $.097\pm.002$ | $-.222\pm.019$ | $.090\pm.005$ | $-.099\pm.054$ | $.095\pm.004$ | $-.112\pm.042$ | $.095\pm.005$ | $-.564\pm.063$
Table 2: Best fit parameters for a linear fit of $\left\langle
g_{c}(\epsilon)\right\rangle$ vs
$\left\langle\bar{\mathscr{A}}(\epsilon)\right\rangle$, across the spectrum
for various interaction strengths. The linear fits, and the uncertainties on
the fit parameters, are computed analogously as in Table 1. Comparing with
Table 1, we see the fit parameters deviate significantly - in particular, the
slope here no longer agrees with 1/12, and the intercepts are closer to zero.

We can derive distributions for $g_{c}(\epsilon)$ at the transition, in
analogy with Section IV. To that end, assume that a generalized relation

$g_{c}(\epsilon)\approx\frac{1}{F}\frac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}\left[-\log\mathcal{A}_{\ell}(\epsilon)-(\ell-1)\alpha\right]$
(39)

holds. As before, let us write
$\mathcal{A}_{\ell}(\epsilon)=\prod_{i=1}^{\ell}\chi_{i}(\epsilon)$, so that:

$g_{c}(\epsilon)\approx\frac{1}{F}\left[-\sum_{i=1}^{N_{\ell}}\frac{N_{\ell}-i+1}{N_{\ell}}\log\chi_{i}(\epsilon)-\frac{N_{\ell}-1}{2}\alpha\right]$

Assuming again $\chi_{i}(\epsilon)\sim\text{Uni}[0,1]$ at the MBL crossover,
$-k\log(\chi_{i}(\epsilon))$ is exponentially distributed with rate parameter
$1/k$, hence $g_{c}$ is (up to the shift $(N_{\ell}-1)\alpha/2F:=\Delta$) a
sum of exponential random variables with rate parameters
$\lambda_{i}=FN_{\ell}/(N_{\ell}-i+1)$ ($i=1,\ldots,N_{\ell}$). Such a sum is
hypoexponentially distributed; the cumulative density for $g_{c}$ is then Cox
(1967):

$P_{F}^{N_{\ell}}(g_{c}\leq
x)=1-\sum_{i=1}^{N_{\ell}}\left(\prod_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N_{\ell}}\frac{\lambda_{j}}{\lambda_{j}-\lambda_{i}}\right)e^{-\lambda_{i}(x+\Delta)}$
(40)

As before, we construct a mixture distribution of the form (18), and fit the
$\beta_{i}$’s by minimizing the squared error (19). Examples at the center of
the spectrum are shown in Fig. 20 and at the band edges in Fig. 21. The
resulting probability densities are shown in Figs. 22 and 23. We note that,
for lower interaction strengths, we find these fits yield comparable or
slightly better (depending on the energy density) squared errors than those of
section IV, whereas for stronger interactions the fits described in the main
text are clearly superior.

(a)

(b)

(c)

(d)

(e)

Figure 20: Numerically observed cumulative density functions of $g_{c}$ for
$\epsilon=0.5$ at the critical point, for interaction strengths a) $U=1$
($W=4.8$) , b) $U=2$ ($W=5.4$), c) $U=3$ ($W=5.9$), d) $U=4$ ($W=6.1$), e)
$U=5$ ($W=6.7$). The blue solid lines overlaid are best fit mixture
distributions of the form (18) (now with $P_{F}^{N_{\ell}}$ now given by
(40)). We see reasonably good agreement, roughly on par or better than of Fig.
10, for $U\leq 3$, and significant deviation for $U>3$. This suggests the
mixture distribution derived from (40) is a good descriptor at low interaction
strengths, whereas that derived from (17) in the main text is far better at
larger interaction strengths. The parameters extracted from the fit are a)
$\beta_{1}=1.0$ b) $\beta_{1}=1.0$, c) $\beta_{1}=0.211$, $\beta_{2}=0.765$,
$\beta_{11}=0.023$ d) $\beta_{1}=1.0$, e) $\beta_{1}=1.0$, and $g_{0}=0$ in
all cases. All quantities are computed in chains of size $L=12$ in the half-
filling sector, and from 10000 disorder realizations. When computing the
squared error (19), we sample $n=100$ points uniformly spaced between the
smallest and largest observed values of $g_{c}$.

(a)

(b)

(c)

(d)

(e)

Figure 21: Numerically observed cumulative density functions of $g_{c}$ for
$\epsilon=0.5$ at the critical point, for interaction strengths a) $U=1$,
$\epsilon=0.8$ ($W=3.0$), b) $U=2$, $\epsilon=0.3$ ($W=4.2$), c) $U=3$,
$\epsilon=0.35$ ($W=5.6$), d) $U=4$, $\epsilon=0.3$ ($W=4.9$), e) $U=5$,
$\epsilon=0.3$ ($W=5.7$). As in Fig. 20, the mixture distributions are
constructed with the $P_{F}^{N_{\ell}}$ given by (40). The energy densities
shown are the energy whose fitted distribution minimized the squared error
(19), for each interaction strength. We can see excellent agreement between
observation and fit at weaker interactions strengths, suggesting that
averaging over the differences $\mathcal{A}_{\ell}(\epsilon)+(\ell-1)\log 2/2$
is an effective descriptor of the physics at the band edges in the weaker
interacting regime. The parameters extracted from the fit are a)
$\beta_{5}=0.603$, $\beta_{6}=0.396$, $\beta_{11}=.001$, b) $\beta_{1}=0.308$,
$\beta_{2}=0.405$, $\beta_{7}=0.189$, $\beta_{8}=0.097$, c) $\beta_{1}=0.095$,
$\beta_{2}=0.056$, $\beta_{3}=0.676$, $\beta_{7}=0.173$, d) $\beta_{1}=0.957$,
$\beta_{11}=0.043$, e) $\beta_{1}=1.0$, and $g_{0}=0$ in all cases. All
quantities are computed exactly as in Fig. 20

(a)

(b)

(c)

(d)

(e)

Figure 22: Numerically observed probability density functions of $g_{c}$ for
$\epsilon=0.5$ at the critical disorder strength, for interaction strengths a)
$U=1$ ($W=4.8$) , b) $U=2$ ($W=5.4$), c) $U=3$ ($W=5.9$), d) $U=4$ ($W=6.1$),
e) $U=5$ ($W=6.7$). The blue solid lines overlaid are the pdfs obtained from
differentiating the best fit mixture distributions shown in Fig. 20.

(a)

(b)

(c)

(d)

(e)

Figure 23: Numerically observed probability density functions of $g_{c}$ at
the critical disorder strength for a) $U=1$, $\epsilon=0.3$ ($W=3.4$), b)
$U=2$, $\epsilon=0.3$ ($W=4.2$), c) $U=3$, $\epsilon=0.35$ ($W=5.6$), d)
$U=4$, $\epsilon=0.3$ ($W=4.9$), e) $U=5$, $\epsilon=0.3$ ($W=5.7$). The blue
solid lines overlaid are the pdfs obtained from differentiating the best fit
mixture distributions shown in Fig. 21.

## References

  * Heußen et al. (2021) S. Heußen, C. D. White, and G. Refael, Phys. Rev. B 103, 064201 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.064201.
  * Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991), URL https://link.aps.org/doi/10.1103/PhysRevA.43.2046.
  * Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994), URL https://link.aps.org/doi/10.1103/PhysRevE.50.888.
  * Rigol et al. (2008) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Rigol and Srednicki (2012) M. Rigol and M. Srednicki, Phys. Rev. Lett. 108, 110601 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.108.110601.
  * Mori et al. (2018) T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 112001 (2018), URL https://doi.org/10.1088/1361-6455/aabcdf.
  * Ueda (2020) M. Ueda, Nature Reviews Physics 2, 668 (2020).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015), eprint https://doi.org/10.1146/annurev-conmatphys-031214-014726, URL https://doi.org/10.1146/annurev-conmatphys-031214-014726.
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958), URL https://link.aps.org/doi/10.1103/PhysRev.109.1492.
  * Huse et al. (2013) D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Phys. Rev. B 88, 014206 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.014206.
  * Pekker et al. (2014) D. Pekker, G. Refael, E. Altman, E. Demler, and V. Oganesyan, Phys. Rev. X 4, 011052 (2014), URL https://link.aps.org/doi/10.1103/PhysRevX.4.011052.
  * Chandran et al. (2014) A. Chandran, V. Khemani, C. R. Laumann, and S. L. Sondhi, Phys. Rev. B 89, 144201 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.144201.
  * Vasseur et al. (2015) R. Vasseur, A. C. Potter, and S. A. Parameswaran, Phys. Rev. Lett. 114, 217201 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.217201.
  * Potter and Vishwanath (2015) A. C. Potter and A. Vishwanath, _Protection of topological order by symmetry and many-body localization_ (2015), URL https://arxiv.org/abs/1506.00592.
  * Bahri et al. (2015) Y. Bahri, R. Vosk, E. Altman, and A. Vishwanath, Nat Commun 6 (2015), URL https://www.nature.com/articles/ncomms8341.
  * Vasseur et al. (2016) R. Vasseur, A. J. Friedman, S. A. Parameswaran, and A. C. Potter, Phys. Rev. B 93, 134207 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.134207.
  * Kells et al. (2018) G. Kells, N. Moran, and D. Meidan, Phys. Rev. B 97, 085425 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.085425.
  * Iadecola and Schecter (2018) T. Iadecola and M. Schecter, Phys. Rev. B 98, 144204 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.144204.
  * Friedman et al. (2018) A. J. Friedman, R. Vasseur, A. C. Potter, and S. A. Parameswaran, Phys. Rev. B 98, 064203 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.064203.
  * Kuno (2019) Y. Kuno, Phys. Rev. Research 1, 032026 (2019), URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.032026.
  * Parameswaran and Vasseur (2018) S. A. Parameswaran and R. Vasseur, Reports on Progress in Physics 81, 082501 (2018), URL https://doi.org/10.1088/1361-6633/aac9ed.
  * Santos et al. (2005) L. F. Santos, M. I. Dykman, M. Shapiro, and F. M. Izrailev, Phys. Rev. A 71, 012317 (2005), URL https://link.aps.org/doi/10.1103/PhysRevA.71.012317.
  * Altshuler et al. (2010) B. Altshuler, H. Krovi, and J. Roland, Proceedings of the National Academy of Sciences 107, 12446 (2010), eprint https://www.pnas.org/doi/pdf/10.1073/pnas.1002116107, URL https://www.pnas.org/doi/abs/10.1073/pnas.1002116107.
  * Laumann et al. (2015) C. Laumann, R. Moessner, A. Scardicchio, and S. Sondhi, The European Physcial Journal Special Topics 224, 75 (2015), URL https://link.springer.com/article/10.1140/epjst/e2015-02344-2#citeas.
  * Yao et al. (2015) N. Y. Yao, C. R. Laumann, and A. Vishwanath, _Many-body localization protected quantum state transfer_ (2015), URL https://arxiv.org/abs/1508.06995.
  * Khemani et al. (2015) V. Khemani, R. Nandkishore, and S. Sondhi, Nat. Phys. 11, 560 (2015), URL https://www.nature.com/articles/nphys3344.
  * Goihl et al. (2020) M. Goihl, N. Walk, J. Eisert, and N. Tarantino, Phys. Rev. Research 2, 013120 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013120.
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.82.174411.
  * Bauer and Nayak (2013) B. Bauer and C. Nayak, 2013, P09005 (2013), URL https://doi.org/10.1088/1742-5468/2013/09/p09005.
  * Žnidarič et al. (2008) M. Žnidarič, T. c. v. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.064426.
  * Kjäll et al. (2014) J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.107204.
  * Luitz et al. (2015) D. J. Luitz, N. Laflorencie, and F. Alet, PhysRevB 91 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.081103.
  * Imbrie (2016) J. Z. Imbrie, J Stat Phys 163 (2016), URL https://link.springer.com/article/10.1007/s10955-01601508-x.
  * Šuntajs et al. (2020) J. Šuntajs, J. Bonča, T. c. v. Prosen, and L. Vidmar, Phys. Rev. E 102, 062144 (2020), URL https://link.aps.org/doi/10.1103/PhysRevE.102.062144.
  * Sels and Polkovnikov (2021) D. Sels and A. Polkovnikov, Phys. Rev. E 104, 054105 (2021), URL https://link.aps.org/doi/10.1103/PhysRevE.104.054105.
  * Abanin et al. (2021) D. Abanin, J. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. Parameswaran, F. Pollmann, A. Potter, M. Serbyn, and R. Vasseur, Annals of Physics 427, 168415 (2021), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/S000349162100021X.
  * Crowley and Chandran (2022) P. J. D. Crowley and A. Chandran, SciPost Phys. 12, 201 (2022), URL https://scipost.org/10.21468/SciPostPhys.12.6.201.
  * Morningstar et al. (2022) A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Phys. Rev. B 105, 174205 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.174205.
  * Huse et al. (2014) D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.174202.
  * Ros et al. (2015) V. Ros, M. Müller, and A. Scardicchio, Nuclear Physics B 891, 420 (2015), ISSN 0550-3213, URL https://www.sciencedirect.com/science/article/pii/S0550321314003836.
  * Serbyn et al. (2013) M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.127201.
  * Chandran et al. (2015) A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin, Phys. Rev. B 91, 085425 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.085425.
  * Rademaker and Ortuño (2016) L. Rademaker and M. Ortuño, Phys. Rev. Lett. 116, 010404 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.010404.
  * Pekker et al. (2017) D. Pekker, B. K. Clark, V. Oganesyan, and G. Refael, Phys. Rev. Lett. 119, 075701 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.075701.
  * Wegner (1994) F. Wegner, Ann. Phys. (Berlin) 506 (1994), URL https://onlinelibrary.wiley.com/doi/10.1002/andp.19945060203.
  * Kehrein (2006) S. Kehrein, _The Flow Equation Approach to Many-Particle Systems_ (Springer Berlin, 2006), ISBN 978-3-540-34067-6.
  * Quito et al. (2016) V. L. Quito, P. Titum, D. Pekker, and G. Refael, Phys. Rev. B 94, 104202 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.104202.
  * Peng et al. (2019) P. Peng, Z. Li, H. Yan, K. X. Wei, and P. Cappellaro, Phys. Rev. B 100, 214203 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.214203.
  * Kulshreshtha et al. (2018) A. K. Kulshreshtha, A. Pal, T. B. Wahl, and S. H. Simon, Phys. Rev. B 98, 184201 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.184201.
  * Varma et al. (2019) V. K. Varma, A. Raj, S. Gopalakrishnan, V. Oganesyan, and D. Pekker, Phys. Rev. B 100, 115136 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.115136.
  * De Roeck and Huveneers (2017) W. De Roeck and F. m. c. Huveneers, Phys. Rev. B 95, 155129 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.155129.
  * Luitz et al. (2017) D. J. Luitz, F. m. c. Huveneers, and W. De Roeck, Phys. Rev. Lett. 119, 150602 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.150602.
  * Thiery et al. (2018) T. Thiery, F. m. c. Huveneers, M. Müller, and W. De Roeck, Phys. Rev. Lett. 121, 140601 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.140601.
  * Thiery et al. (2017) T. Thiery, M. Müller, and W. D. Roeck (2017), eprint 1711.09880.
  * Johri et al. (2015) S. Johri, R. Nandkishore, and R. N. Bhatt, Phys. Rev. Lett. 114, 117401 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.117401.
  * Panda and Banerjee (2020) A. Panda and S. Banerjee, Phys. Rev. B 101, 184201 (2020), URL https://link.aps.org/doi/10.1103/PhysRevB.101.184201.
  * Hatano and Nelson (1996) N. Hatano and D. R. Nelson, Phys. Rev. Lett. 77, 570 (1996), URL https://link.aps.org/doi/10.1103/PhysRevLett.77.570.
  * Hatano and Nelson (1998) N. Hatano and D. R. Nelson, Phys. Rev. B 58, 8384 (1998), URL https://link.aps.org/doi/10.1103/PhysRevB.58.8384.
  * Hatano and Nelson (1997) N. Hatano and D. R. Nelson, Phys. Rev. B 56, 8651 (1997), URL https://link.aps.org/doi/10.1103/PhysRevB.56.8651.
  * Kato (1995) T. Kato, _Perturbation theory in a finite-dimensional space_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1995), pp. 62–126, ISBN 978-3-642-66282-9, URL https://doi.org/10.1007/978-3-642-66282-9_2.
  * Heiss and Sannino (1990) W. D. Heiss and A. L. Sannino, Journal of Physics A: Mathematical and General 23, 1167 (1990), URL https://doi.org/10.1088%2F0305-4470%2F23%2F7%2F022.
  * Heiss and Sannino (1991) W. D. Heiss and A. L. Sannino, Phys. Rev. A 43, 4159 (1991), URL https://link.aps.org/doi/10.1103/PhysRevA.43.4159.
  * Heiss (2000) W. D. Heiss, Phys. Rev. E 61, 929 (2000), URL https://link.aps.org/doi/10.1103/PhysRevE.61.929.
  * Heiss (2001) W. D. Heiss, AIP Conference Proceedings 597, 311 (2001), eprint https://aip.scitation.org/doi/pdf/10.1063/1.1427477, URL https://aip.scitation.org/doi/abs/10.1063/1.1427477.
  * Heiss and Harney (2001) W. D. Heiss and H. Harney, Eur. Phys. J. D 17, 149 (2001), URL https://link.springer.com/article/10.1007/s100530170017.
  * Heiss (2004) W. D. Heiss, Journal of Physics A: Mathematical and General 37, 2455 (2004), URL https://doi.org/10.1088/0305-4470/37/6/034.
  * Heiss (2008) W. D. Heiss, Journal of Physics A: Mathematical and Theoretical 41, 244010 (2008), URL https://doi.org/10.1088/1751-8113/41/24/244010.
  * Heiss (2012) W. D. Heiss, Journal of Physics A: Mathematical and Theoretical 45, 444016 (2012), URL https://doi.org/10.1088/1751-8113/45/44/444016.
  * Luitz and Piazza (2019) D. J. Luitz and F. Piazza, Phys. Rev. Research 1, 033051 (2019), URL https://link.aps.org/doi/10.1103/PhysRevResearch.1.033051.
  * Hamazaki et al. (2019) R. Hamazaki, K. Kawabata, and M. Ueda, Phys. Rev. Lett. 123, 090603 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.090603.
  * Feinberg and Zee (1999a) J. Feinberg and A. Zee, Phys. Rev. E 59, 6433 (1999a), URL https://link.aps.org/doi/10.1103/PhysRevE.59.6433.
  * Feinberg and Zee (1999b) J. Feinberg and A. Zee, Nuclear Physics B 552, 599 (1999b), ISSN 0550-3213, URL https://www.sciencedirect.com/science/article/pii/S0550321399002461.
  * Brézin and Zee (1998) E. Brézin and A. Zee, Nuclear Physics B 509, 599 (1998), ISSN 0550-3213, URL https://www.sciencedirect.com/science/article/pii/S0550321397006524.
  * Potirniche et al. (2019) I.-D. Potirniche, S. Banerjee, and E. Altman, Phys. Rev. B 99, 205149 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.205149.
  * Burnham and Anderson (2004) K. P. Burnham and D. R. Anderson, Sociological Methods & Research 33, 261 (2004), eprint https://doi.org/10.1177/0049124104268644, URL https://doi.org/10.1177/0049124104268644.
  * Villalonga and Clark (2020a) B. Villalonga and B. K. Clark, _Characterizing the many-body localization transition through correlations_ (2020a), URL https://arxiv.org/abs/2007.06586.
  * Villalonga and Clark (2020b) B. Villalonga and B. K. Clark, _Eigenstates hybridize on all length scales at the many-body localization transition_ (2020b), URL https://arxiv.org/abs/2005.13558.
  * Rodriguez et al. (2011) A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. Römer, Phys. Rev. B 84, 134209 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.134209.
  * Cox (1967) D. R. Cox, _Renewal theory._ , Methuen’s monographs on applied probability and statistics (Methuen; [distributed in the USA by Barnes & Noble, 1967), URL https://search.ebscohost.com/login.aspx?direct=true&AuthType=ip,sso&db=cat09073a&AN=cit.oai.edge.caltech.folio.ebsco.com.fs00001057.62fa39eb.c0bc.4443.8272.a7be9a56a9fe&site=eds-live&scope=site.