# Localization of a mobile impurity interacting with an Anderson insulator

Pietro Brighi1 Alexios A. Michailidis1,2 Kristina Kirova1 Dmitry A. Abanin2
Maksym Serbyn1 1IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria
2Department of Theoretical Physics, University of Geneva, 24 quai Ernest-
Ansermet, 1211 Geneva, Switzerland

(August 27, 2024)

###### Abstract

Thermalizing and localized many-body quantum systems present two distinct
dynamical phases of matter. Recently, the fate of a localized system coupled
to a thermalizing system viewed as a quantum bath received significant
theoretical and experimental attention. In this work, we study a mobile
impurity, representing a small quantum bath, that interacts locally with an
Anderson insulator with a finite density of localized particles. Using static
Hartree approximation to obtain an effective disorder strength, we formulate
an analytic criterion for the perturbative stability of the localization.
Next, we use an approximate dynamical Hartree method and the quasi-exact time-
evolved block decimation (TEBD) algorithm to study the dynamics of the system.
We find that the dynamical Hartree approach which completely ignores
entanglement between the impurity and localized particles predicts the
delocalization of the system. In contrast, the full numerical simulation of
the unitary dynamics with TEBD suggests the stability of localization on
numerically accessible timescales. Finally, using an extension of the density
matrix renormalization group algorithm to excited states (DMRG-X), we
approximate the highly excited eigenstates of the system. We find that the
impurity remains localized in the eigenstates and entanglement is enhanced in
a finite region around the position of the impurity, confirming the dynamical
predictions. Dynamics and the DMRG-X results provide compelling evidence for
the stability of localization.

##  I Introduction

Many body localization (MBL) is an example of a stable dynamical phase of
matter that avoids thermal equilibrium. The existence of the MBL phase may be
understood from the stability of localization in the Anderson Hamiltonian [1]
with respect to weak but finite interactions [2, 3, 4]. Recently a number of
dynamical properties of MBL [5, 6] were explained via the existence of an
extensive number of quasilocal conserved quantities [7, 8, 9]. At the same
time, experiments provided strong support for the stability of the MBL phase
on long timescales, verified a number of theoretical predictions, and started
exploring new regimes where the theoretical understanding is incomplete [10,
11, 12, 13]. In particular, experiments suggested the existence of MBL in
higher dimension [14, 15] and also probe the so-called many-body mobility
edges [16], although theoretically their existence is subject to debate [17,
18, 19, 20, 21, 22, 23, 24, 25, 26].

Another open question concerns the stability of localization in an MBL system
coupled to a so-called quantum bath, represented by another quantum system
that thermalizes in the absence of the coupling. In the case of a
thermodynamically large bath, where the back-action of the localized system
can be neglected, the system-bath coupling is expected to result in
delocalization. However, considering a bath whose dimension is comparable to
the localized system, or smaller, can yield distinct outcomes, especially if
the back-action on the bath is taken into account. The MBL degrees of freedom
can localize the bath – a phenomenon dubbed “MBL proximity effect” by
Nandkishore [27]. Alternatively, the bath can thermalize the formerly
localized system [28].

Inspired by the avalanche mechanism [22] for the delocalization transition, a
number of works studied the effect of a “thermal grain” coupled to a localized
system [29, 28, 30]. These works described the bath by a random matrix type
Hamiltonian and account for the quasi-local nature of integrals of motion
while considering coupling between the bath and localized system. In these
models the bath lacks spatial structure, thus it cannot be localized by the
MBL system, excluding _a priori_ the MBL proximity effect. In order to keep
the microscopic structure of the bath, Refs. [27, 31, 32] represented it by a
set of thermalizing particles. In this framework, thermal and localized
degrees of freedom coexist and are coupled through local interaction. It has
been numerically shown [31] that localization can globally persist, if the
bandwidth of the thermal particles is small. However, the fingerprints of
localization of the bath under the influence of the disordered degrees of
freedom was not studied in detail and still remains an open question.

All the studies discussed so far relied on the use of exact diagonalization
(ED), which dramatically limits the system sizes available. In contrast,
recent experimental studies [33, 34] have addressed this problem on bosonic
quantum simulators, enabling the study of large systems. In Ref. [33] a
“global” setup was used, where the thermal degrees of freedom are
homogeneously distributed through the $2$d lattice. Varying the number of
thermal particles, the experiment showed evidence of the stability of
localization when the thermalizing bosons are a small fraction of the total.
On the other hand, the authors of Ref. [34] used a different approach. There,
a $1$d chain is split into a disorder-free segment, that represents a bath,
and is connected to a disordered segment. The experiment investigated the
stability of localization while changing the size of the disorder-free
segment. While localization is stable for small thermal chains, signs of
delocalization were observed when the bath constitutes half of the whole
system.

In this work, inspired by Ref. [33], we study a one-dimensional system of two
hard-core bosonic species interacting with one another. The disordered
particles form an Anderson insulator and are coupled locally to a single
disorder-free boson representing the smallest possible bath. Analyzing the
model under a mean-field type approximation, we formulate a perturbative
criterion for stability of localization, that suggests that localization
remains robust at strong interaction and disorder. This is in contrast with
the results discussed in Ref. [35], where ergodicity is introduced in the
system by a stable doublon, that can freely move in the strong interaction
limit, therefore facilitating delocalization in the case when single-particle
localization length is long enough.

In an accompanying paper [36], we studied the dynamics of the model discussed
here in large systems using matrix product state (MPS) [37] techniques. There,
at strong interactions and disorder we observed the localization of the bath
through the back-action of the Anderson insulator. Here, besides analytic
results supporting localization in such a regime, we provide a thorough
comparison of fully interacting dynamics with an approximate time-dependent
Hartree method, which neglects entanglement and quantum correlations among the
two particle species. This latter technique reveals delocalization at long
times, supported by the diffusive spreading of the particle constituting the
bath and decaying memory of the initial state of the Anderson insulator. This
result highlights the fundamental role played by entanglement, which we study
in detail in the present work.

Finally, we use density-matrix renormalization group for excited states method
(DMRG-X) [38, 39, 40], to study highly excited eigenstates of the Hamiltonian
in large chains. This effectively allows us to probe localization at infinite
times. Analyzing the expectation value of density in eigenstates, we observe
localization of the small bath due to the interaction with the Anderson
insulator. Furthermore, we find that eigenstates show area-law entanglement,
thus providing complementary support for the persistence of localization at
strong interactions.

The structure of this paper is as follows. In Section II we introduce the
model and describe the typical initial state chosen for the dynamics. In
Section III we analyze the Hartree limit and study perturbatively the effect
of the interaction. Following that, in Section IV, we introduce the time-
dependent Hartree approximation for the dynamics, showing its results and
comparing them with quasi-exact TEBD results. Finally, in Section V we present
our DMRG-X study.

##  II Model

In this work, we study two bosonic species of particles, interacting through
an on-site potential. The disordered bosons ($d$-bosons) are subject to a
random potential drawn from a uniform distribution $\epsilon_{i}\in[-W,W]$ and
are governed by the Hamiltonian $\hat{H}_{d}$

$\hat{H}_{d}=t_{d}\sum_{i=1}^{L-1}(\hat{d}^{\dagger}_{i}\hat{d}_{i+1}+\text{h.c.})+\sum_{i=1}^{L}\epsilon_{i}\hat{n}_{d,i},$
(1)

where $\hat{d_{i}}$ is the annihilation operator for the $d$-bosons,
$\hat{n}_{d,i}=\hat{d}^{\dagger}_{i}\hat{d}_{i}$ is their density operator,
and $t_{d}$ is the hopping strength. These bosons realize an Anderson
insulator [1], that is localized at any density of bosons thereby providing a
specific system that avoids thermalization.

The small quantum bath will be represented by the clean bosons ($c$-bosons)
described by $\hat{H}_{c}$,

$\hat{H}_{c}=t_{c}\sum_{i=1}^{L-1}(\hat{c}^{\dagger}_{i}\hat{c}_{i+1}+\text{h.c.}),$
(2)

that are characterized by a single hopping parameter $t_{c}$ and are not
subject to disordered potential. Finally, the two boson species are coupled
via the on-site Hubbard interaction,

$\displaystyle\hat{H}_{\text{int}}=U\sum_{i=1}^{L}\hat{n}_{c,i}\hat{n}_{d,i},$
(3)

where $U$ is the interaction strength and
$\hat{n}_{c,i}=\hat{c}^{\dagger}_{i}\hat{c}_{i}$ is the number operator of
$c$-bosons.

The full interacting Hamiltonian then reads:

$\hat{H}=\hat{H}_{d}+\hat{H}_{c}+\hat{H}_{\text{int}}.$ (4)

The system described by the Hamiltonian (4) has $U(1)\times U(1)$ symmetry, as
the number of both types of particles,
$\hat{N}_{c/d}=\sum_{i}\hat{n}_{c/d,i}$, is conserved. In what follows we
restrict to a sector of finite $d$-bosons density and to the case
$t_{d}=t_{c}=1$. In particular, density of $\nu_{d}=N_{d}/L=1/3$ would be used
unless specified otherwise. At the same time, throughout this work we consider
the presence of a _single_ $c$-boson, that realizes the smallest possible
quantum bath with non-trivial spatial structure and local coupling to a
localized system. In this respect, our setting resembles the experimental
setup of Ref. [33] which was, however, performed on a two dimensional lattice
and considered various densities of clean particles.

Although the boson density is sufficient to specify a particular sector of the
Hilbert space, we further restrict ourselves to states where $d$-bosons have a
globally homogeneous distribution. In particular, we will study initial states
corresponding to a $d$-bosons density wave, with the single $c$-boson located
on the central site, as exemplified by $\ket{\psi_{0}}$ below on a system of
$L=18$ sites and with $\nu_{d}=1/3$

$\ket{\psi_{0}}=|{\bullet}{\circ}{\circ}{\bullet}{\circ}{\circ}{\bullet}{\circ}{\color[rgb]{1,0,0}\bullet}{\bullet}{\circ}{\circ}{\bullet}{\circ}{\circ}{\bullet}{\circ}{\circ}\rangle,$
(5)

where empty circles represent empty sites, and black (red) circles are sites
occupied by $d$\- and $c$-boson respectively. Such initial state resembles the
configurations used in experiments and can be characterized by the so-called
imbalance [11, 41], quantifying the memory of the initial density-wave
configuration in the system after a quench.

While a strictly periodic arrangement of $d$-bosons akin to the state (5) is
not required, we assume that the density of $d$-bosons is on average
distributed uniformly on a scale that is larger than a few lattice spacings.
This assumption is important since the presence of large empty/occupied
regions in the chain would imply the effective absence of disorder for
$c$-boson in that region. While such configurations could be used to imitate
another experimental study of MBL-bath coupling [34], states where extensive
regions of the chain are fully empty or occupied by $d$-bosons are far from
typical initial product states. Moreover, we expect that if the localization
of $d$-bosons persists, such states are nearly decoupled from spatially
homogenous initial states of the type (5). Indeed, in order to connect the
highly inhomogeneous state such as
$\ket{\psi_{0}}=|{\bullet}{\bullet}{\bullet}{\bullet}{\bullet}{\bullet}{\circ}{\circ}{\circ}{\circ}{\circ}{\color[rgb]{1,0,0}\bullet}{\circ}{\circ}{\circ}{\circ}{\circ}{\circ}\rangle$
to the density-wave state in Eq. (5), the tunneling of an extensive number of
$d$-boson over long distances is required.

Figure 1:  (a) Subtracting deterministic components from the distribution of
$\langle\hat{n}_{d,i}\rangle$ (inset) results in an approximately Gaussian
distribution of the effective disorder potential $\tilde{\epsilon}_{i}$, shown
here for $\nu_{d}=1/3$. Standard deviation of this distribution, $\sigma$, is
used to define the effective disorder strength, $\tilde{W}/U=\sigma$. (b) The
effective disorder strength $\tilde{W}$ shows a non-monotonic behavior as a
function of $W$ with a maximum at $W^{*}\approx 5$. The data shown are
obtained for $L=500$ sites and averaged over $50$ disorder realizations.

##  III Hartree approximation and criterion for stability of localization

First, we formulate the Hartree approximation and use it to study the
Hamiltonian (4). This facilitates the choice of parameters in the Hamiltonian,
and allows to formulate an analytic criterion for stability of localization
with respect to two-particle tunneling processes.

###  III.1 Effective Disorder

The Hartree approximation adopted in this section consists of replacing the
density operator $\hat{n}_{d,i}$ in Eq. (4) with its infinite-time average
$\langle\hat{n}_{d,i}\rangle$. This approximation would be fully justified in
the case of instantaneous relaxation of $d$-bosons with respect to $c$-boson
dynamics, i.e. if $t_{d},W\gg t_{c}$. Moreover, as in an Anderson insulator
the fluctuations of the expectation value $\langle\hat{n}_{d,i}\rangle$ are
finite at all times, the Hartree approximation overestimates localization. In
spite of these shortcomings, the Hartree approximation will assist us with the
choice of model parameters and will also allow defining an effective disorder
strength thus quantifying the analytic criterion for stability of
localization.

The infinite-time average of the $d$-bosons density is given by the diagonal
ensemble generated by the eigenstates of $\hat{H}_{d}$:

$\begin{split}\langle\hat{n}_{d,i}\rangle&=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}dt\,\bra{\psi(t)}\hat{n}_{d,i}\ket{\psi(t)}\\\
&=\sum_{n}|c_{n}|^{2}\bra{E_{n}}\hat{n}_{d,i}\ket{E_{n}},\end{split}$ (6)

where $\\{\ket{E_{m}}\\}$ is the eigenbasis of $\hat{H}_{d}$,
$c_{n}=\langle\psi_{0}\ket{E_{n}}$ and
$\ket{\psi(t)}=\hat{T}e^{-\imath\int_{0}^{t}\hat{H}_{d}(t^{\prime})dt^{\prime}}\ket{\psi_{0}}=e^{-\imath\hat{H}_{d}t}\ket{\psi_{0}}$.
The initial state $\ket{\psi_{0}}$ is taken to be a density wave, see Eq. (5).
Due to the random potential in the Anderson Hamiltonian, the $d$-bosons
occupation at infinite times acquires a quasi-random nature and is distributed
according to $P(\langle\hat{n}_{d,i})\rangle$ in the range $[0,1]$.

Thanks to the noninteracting nature of the problem in the Hartree
approximation, the eigenstates in Eq. (6) can be written as product states of
single particle orbitals $\ket{\ell}$. Then, by going into eigenbasis one can
decouple the sum over ${L}\choose{N_{d}}$ terms in Eq. (6) into a double sum
over the $N_{d}$ occupied sites in the initial state and the $L$ single
particle orbitals. The problem of finding the infinite time occupation values
then reduces to finding eigenstates of a single-particle Hamiltonian, whose
Hilbert space scales as $L$, thus greatly reducing the complexity of the
problem. The infinite-time average of the $d$-bosons density, hence, can be
carried out through Eq. (6) for very large systems, $L=500$, where boundary
effects on its probability distribution are negligible, and at various
disorder strengths.

The infinite-time occupation of the $d$-bosons plays the role of an effective
random chemical potential experienced by the $c$-boson in the Hartree
approximation. As shown in the inset of Figure 1(a), when the disorder is
strong compared to the hopping parameter, $t_{d}=1$,
$P(\langle\hat{n}_{d,i})\rangle$ resembles a bimodal distribution, with two
distinct peaks at $\langle\hat{n}_{d,i}\rangle\sim 0,1$. This behavior can be
understood as a consequence of the strong localization of $d$-bosons for
disorder $W=10$. Thus, even at infinite time, the $d$-bosons remain close to
their initial position and expectation value of $\langle\hat{n}_{d,i}\rangle$
remains close to its initial value, i.e. either $0$ or $1$, depending on the
considered site. As disorder is decreased, the two peak structure gradually
disappears, and the distribution acquires a single peak around the average
density $\nu_{d}$.

In order to define the effective disorder $\tilde{\epsilon}_{i}$ generated by
the $d$-bosons, we need to isolate the random part of the distribution of
$\langle\hat{n}_{d,i}\rangle$. To this end, we subtract from
$\langle\hat{n}_{d,i}\rangle$ the uniform and the period-$1/\nu_{d}$
contribution. In Fourier space this corresponds to modifying the Fourier
harmonics of density, $\tilde{n}_{d}(k)=\sum_{j}\langle\hat{n}_{d,j}\rangle
e^{-\imath kj}$, as follows:

$\tilde{\epsilon}(k)=\tilde{n}_{d}(k)-(1-\alpha)L\nu_{d}\delta_{k,0}-\alpha\sum_{j=1}^{L}f(j,\nu_{d})e^{-\imath
kj},$ (7)

where $\alpha=\sum_{j}\langle\hat{n}_{d,j}\rangle e^{-\imath\pi\nu_{d}j}$
corresponds to the weight of the $1/\nu_{d}$ harmonic in the particle
distribution and $f(j,\nu_{d})=\sum_{n=0}^{N_{d}-1}\delta_{j,n/\nu_{d}+1}$ is
a functional representation of the initial density wave configuration.
Transforming $\tilde{\epsilon}(k)$ back to the real space, we obtain the
effective random potential $\tilde{\epsilon}_{i}$. Its distribution differs
from the one of the infinite-time density, as it can be seen in Figure 1(a).
In particular, $P(\tilde{\epsilon}_{i})$ is centered around zero for all
values of disorder, and has an approximately Gaussian shape. Thus, we use the
standard deviation of this distribution as the effective disorder strength
$\tilde{W}=U\mathop{\rm std}P(\tilde{\epsilon}_{i})$ experienced by $c$-boson
in the Hartree approximation.

The effective disorder strength measured in units of $U$, $\tilde{W}/U$ is
shown as a function of the disorder strength experienced by $d$-bosons, $W$,
in Figure 1(b). First, we note a non-monotonic dependence of $\tilde{W}$ on
the $d$-boson disorder strength, $W$. The effective disorder strength
$\tilde{W}$ presents a maximum around $W^{*}\approx 5$, whose position depends
weakly on the density of $d$-bosons. This behavior can be naturally explained
by considering two opposite limits: at weak $W$ the localization length of
$d$-bosons is much larger and the initial period-2/3/4 density wave
configuration is washed out at late times resulting in weak effective disorder
$\tilde{W}$. In the opposite limit of very strong $W$, the $d$-bosons remain
frozen close to their initial positions resulting in a nearly perfect periodic
potential experienced by $c$-boson. However such periodic potential is unable
to localize $c$-boson and is subtracted in Eq. (7), thus again resulting in a
weak effective disorder. Given that $\tilde{W}$ is expected to decrease for
very large and small $W$, we expect it to achieve a maximal value at some
intermediate disorder $W$. Finally, we study in Fig. 1(b) the dependence of
the effective disorder on the $d$-bosons density $\nu_{d}$. As $\nu_{d}$ is
increased, $\tilde{W}$ increases accordingly due to the fact that the
effective random potential $\tilde{\epsilon}_{i}$ is generated by the
$d$-bosons density. We note, however, that for $\nu_{d}>1/2$ the effective
disorder would decrease again, because of the hard-core nature of the bosons.

###  III.2 Localization length of the $c$-boson in Hartree approximation

We demonstrated that in the Hartree approximation the $c$-boson experiences an
approximately Gaussian-distributed random potential with disorder strength
$\tilde{W}$ that depends on the initial state and disorder experienced by
$d$-bosons. Since an arbitrary weak disorder potential suffices to localize a
single particle in a one-dimensional lattice, the $c$-boson in the Hartree
approximation is always localized. We proceed with the calculation of its
localization length $\xi_{c}$ that provides a characteristic lengthscale of
localization and can be compared with the localization length of the
$d$-bosons, $\xi_{d}$.

To obtain the localization length for the $c$-boson we use the weak disorder
approximation, justified at weak values of $U$ as the random potentials
$\tilde{\epsilon}_{i}$ are restricted to $[-1,1]$. Perturbing around the
tight-binding limit with the transfer matrix method [42], we obtain
$\xi_{c}(k)\approx{8t^{2}_{c}\sin^{2}(k)}/{\langle\tilde{\epsilon}^{2}_{i}\rangle}$
that depends on the momentum $k$ that determines single particle energy in
absence of disorder, $E(k)=-2t_{c}\cos(k)$. The average localization length is
calculated by performing the integral over the complete band

$\xi_{c}\approx\frac{8t^{2}_{c}}{\langle\tilde{\epsilon}^{2}_{i}\rangle}\int_{-2t_{c}}^{2t_{c}}d\varepsilon\left[1-\left(\frac{\varepsilon}{2t_{c}}\right)^{2}\right]\rho(\varepsilon),$
(8)

where $\rho(\varepsilon)$ is the usual density of states of the tight-binding
Hamiltonian. Carrying out the integral we realize that it contributes only to
a numerical factor as the $c$-boson hopping terms cancel out. Recalling that,
since $\langle\tilde{\epsilon}_{i}\rangle=0$, the variance
$\langle\tilde{\epsilon}^{2}_{i}\rangle$ corresponds to the definition of the
effective disorder strength squared, we obtain

$\xi_{c}\approx\frac{4t_{c}^{2}}{a\tilde{W}^{2}},$ (9)

where $a=1$ is the lattice spacing.

Figure 2:  The localization length of the $c$-boson from the Hartree
approximation for initial period-3 density wave state of $d$-bosons. Inset
shows bare data, for different values of $U$ ranging from $U=0.5$ (dark blue)
to $U=5$ (yellow). For a broad range of interaction strength, $\xi_{c}$ has a
minimum at $W^{*}\approx 5$. The main plot shows the collapse of data,
confirming the scaling $\xi_{c}\sim U^{-2}$.  Figure 3:  (a) The typical ratio
between matrix element and level spacing as a function of disorder for
different coupling strengths $U$ rapidly decreases as the interaction strength
increases. In the left inset we compare analytic expression for $\mathcal{R}$
with its numerical estimate. The right inset shows the tunneling process
induced by $\hat{H}_{\text{int}}$. (b) Phase diagram resulting from the
criterion ${\cal R}=1$ reveals a broad region where the localization is
perturbatively stable, which increases with increasing $U$ and $W$.

The resulting simple expression for $\xi_{c}$, Eq. (9), has two main
consequences. First, we expect that $\xi_{c}$ inherits the non-monotonic
dependence on $W$ and has a minimum approximately when the effective disorder
strengths is maximal in Fig. 1. Second, since $\tilde{W}\propto U$, we expect
that localization length diverges as $\xi_{c}\propto 1/U^{2}$ at small values
of $U$ where weak disorder approximation is controllable. To confirm these
predictions, we extract the localization length of $c$-bosons from numerical
simulations. Using the effective disorder obtained from the eigenstates of the
$d$-boson Hamiltonian, we use exact diagonalization to calculate the single-
particle wave functions of $c$-bosons resulting from the Hamiltonian
$H_{c}^{\text{Hartree}}=H_{c}+\sum_{i}U\langle\hat{n}_{d,i}\rangle\hat{n}_{c,i}$.
Each eigenstate $\ket{\phi_{\varepsilon}}$ of this Hamiltonian can be
characterized by an (energy-dependent) localization length
$\xi_{c}(\varepsilon)$. After obtaining $\xi_{c}(\varepsilon)$ for each
eigenstate through an exponential fit of its probability distribution in the
real space, $|\bra{i}\phi_{\varepsilon}\rangle|^{2}$, we average the
localization length over all eigenstates and further over $50$ disorder
realizations.

The localization length resulting from the numerical simulation is shown in
Fig. 2, where the non-monotonic behavior of $\xi_{c}$ and its scaling with the
interaction strength becomes apparent. We note that for the adopted choice of
hoppings $t_{d}=t_{c}=1$, $\xi_{c}$ always exceeds the localization length of
$d$-bosons and is tunable by the interaction strength in a broad range. In
particular, for disorder strength $W=6.5$, used in the accompanying paper
[36], $\xi_{d}\approx 0.5$ is smaller than one lattice spacing, while
$\xi_{c}$ takes values from $1.5$ to $100$ lattice spacings, as the
interaction strength is decreased from $U=5$ to $U=0.5$ respectively. This
result also assists in the choice of disorder strength: in order to facilitate
the numerical studies, we fix disorder at $W=6.5$ so that the localization
length of $c$-bosons is close to its minimum.

###  III.3 Analytic criterion for stability of localization

The drawback of the Hartree approximation presented above is that it ignores
fluctuations of $d$-bosons density, thus strongly favoring localization.
However, it provides a useful starting point to address the perturbative
stability of such localized system. To this end, we use the basis of Anderson
localized orbitals provided by Hartree approximation to address the stability
of the system with respect to interactions between clean and disordered
bosons.

We transform the interaction Hamiltonian (3) to the basis of localized
orbitals via the relation
$\hat{d}_{\alpha}=\sum_{i}\psi_{d\alpha}(i)\hat{d}_{i}$, where
$\hat{d}_{\alpha}$ is the annihilation operator of an Anderson orbital and
$\psi_{d\alpha}(i)=\langle{\alpha}|i\rangle$ is the corresponding wave
function and similar expressions hold for $c$-bosons. This leads to

$\hat{H}_{\text{int}}=\sum_{\alpha\beta\gamma\delta}V_{\alpha\beta\gamma\delta}\hat{d}^{\dagger}_{\alpha}\hat{c}^{\dagger}_{\beta}\hat{d}_{\gamma}\hat{c}_{\delta}.$
(10)

As shown in the schematic representation in the inset of Figure 3, the matrix
element $V_{\alpha\beta\gamma\delta}$ corresponds to a correlated hopping
process with its value being given by the overlap of the wave functions of the
corresponding localized orbitals,

$V_{\alpha\beta\gamma\delta}=U\sum_{i=1}^{L}\psi^{*}_{d\alpha}(i)\psi^{*}_{c\beta}(i)\psi_{d\gamma}(i)\psi_{c\delta}(i),$
(11)

where the envelope of the wave functions decays on the scale of the
corresponding localization length, $|\psi_{d\alpha}(i)|\sim
e^{-|i-x_{\alpha}|/(2\xi_{d})}/{\sqrt{2\xi_{d}}}$, where $x_{\alpha}$ is the
site around which the orbital is localized. In order to address the stability
of localization with respect to such correlated hoppings triggered by
interaction, we compare the matrix element of this process to the
corresponding level spacing.

Keeping only the leading behavior of the wave functions, i.e. the exponential
decay, and neglecting their oscillatory behavior results in an upper bound for
the matrix element in Eq. (11) and thus favors delocalization. In this case,
the matrix element can be easily estimated as

$V_{\alpha\beta\gamma\delta}\approx\frac{U}{4\xi_{c}\xi_{d}}\sum_{i=1}^{L}e^{-\frac{|x_{\alpha}-i|+|x_{\gamma}-i|}{\xi_{d}}}e^{-\frac{|x_{\beta}-i|+|x_{\delta}-i|}{\xi_{c}}},$
(12)

from where one can immediately realize that it is exponentially suppressed
whenever any of terms in the exponents $|x-i|/\xi_{c/d}>1$ exceeds one. This
restricts the localized orbitals that can efficiently participate in the
hopping process:

$\begin{cases}|x_{\alpha}-x_{\gamma}|&\lesssim\xi_{d},\quad|x_{\alpha}-x_{\beta}|\lesssim\max(\xi_{c},\xi_{d})\\\
|x_{\beta}-x_{\delta}|&\lesssim\xi_{c},\quad|x_{\gamma}-x_{\delta}|\lesssim\max(\xi_{c},\xi_{d}).\end{cases}$
(13)

Using the fact that $\xi_{c}\geq\xi_{d}$, the summation cancels the $\xi_{d}$
in the denominator resulting in a simple estimate

$V_{\alpha\beta\gamma\delta}\approx\frac{U}{4\xi_{c}}.$ (14)

The second element needed to understand if such tunneling processes are
resonant is the typical level spacing $\delta_{cd}$, obtained as the ratio of
the typical energy difference $\Delta E$ to the number of states $\mathcal{N}$
connected by the interaction. Since the two states differ for the position of
the two bosons, the typical energy difference is given by the sum of the two
disorder strengths: $\Delta E\approx W+\tilde{W}$. To account for the number
of states $(\ket{\psi},\ket{\psi^{\prime}})$ connected by the interaction we
need to consider the constraints (13), which restrict the possible
configurations. First, in the initial state the two bosons must lie within a
distance $\xi_{c}$ from one another, as $\xi_{c}\geq\xi_{d}$, thus
contributing $2\xi_{c}$ possible configurations. Due to the localized nature
of wave functions, the position of the bosons in the final state must lie
within a distance $\xi_{c}$ and $\xi_{d}$ for the $c$-boson and the $d$-bosons
respectively. A naïve computation would then give
$\mathcal{N}=(2\xi_{c})(2\xi_{c})(2\xi_{d})=8\xi^{2}_{c}\xi_{d}$. However, one
must be careful as also in the final state the two bosons must be separated by
at most $\xi_{c}$ for the matrix element not to vanish. This reduces the total
number of states, as not all the moves increasing the distance among the
bosons are allowed. A careful computation reveals that for a given $d$-boson
hopping there are only $(3/2)\xi_{c}$ possible choices, thus reducing the
total number of states. Finally, one has to take into account the finite
density of the $d$-bosons, together with their hard-core nature, that requires
that an initially occupied site must be left empty after the hopping and vice
versa, leading to

$\mathcal{N}\approx
6\xi^{2}_{c}\xi_{d}\nu_{d}(1-\nu_{d}),\quad\delta_{cd}=(W+\tilde{W})/\mathcal{N}.$
(15)

Gathering the results of Eqs. (14)-(15) we obtain the expression for the ratio
of level spacing to the matrix element,

$\mathcal{R}=\frac{3}{2}U\frac{\xi_{c}\xi_{d}\nu_{d}(1-\nu_{d})}{W+\tilde{W}}.$
(16)

Condition $\mathcal{R}<1$ provides a criterion for stability of localization.
Recalling that the $c$-boson localization length depends on the interaction
strength $\xi_{c}\approx\xi_{c}^{0}/U^{2}$, the criterion $\mathcal{R}<1$ can
be rearranged as

$(W+\tilde{W})U>\frac{3}{2}\xi_{c}^{0}\xi_{d}\nu_{d}(1-\nu_{d}),$ (17)

suggesting that localization remains stable at strong interactions and large
disorder.

The typical probability of resonance $\mathcal{R}$ can also be evaluated
numerically and compared with the prediction of Eq. (16). To this end, after
diagonalizing the Hamiltonian in the Hartree approximation, for each $c$-boson
eigenstate we select states lying within the localization lengths $\xi_{c}$
and $\xi_{d}$ from one another and evaluate for all of them the matrix element
$V_{\alpha\beta\gamma\delta}$ and the energy difference. Then we define the
resonance value $\mathit{r}_{\beta}$ fixing the initial state of the
$c$-boson, e.g. $\beta$, and maximizing the ratio of the matrix element to the
energy difference over all other available states. In a similar way we define
the resonance $\mathit{r}_{\alpha\beta}$, where both the initial state of the
$c$-boson ($\beta$) and of the $d$-boson ($\alpha$) are fixed. Being a more
constrained version of $\mathit{r}_{\beta}$, the following inequality holds

$\mathit{r}_{\beta}=\max_{\alpha\gamma\delta}\Bigr{[}\frac{V_{\alpha\beta\gamma\delta}}{\Delta
E}\Bigr{]}\geq\mathit{r}_{\alpha\beta}=\max_{\gamma\delta}\Bigr{[}\frac{V_{\alpha\beta\gamma\delta}}{\Delta
E}\Bigr{]}.$ (18)

Finally, we obtain the typical resonance probability by taking the median of
the distribution of $\mathit{r}_{\beta}$ and $\mathit{r}_{\alpha\beta}$. The
results of this numerical evaluation are shown in the inset of Figure 3(a)
together with the analytic prediction (16) and show an overall agreement. In
particular, we notice that $\mathit{r}_{\alpha\beta}$ is always smaller than
$\mathit{r}_{\beta}$.

From Figure 3(a) we observe that the disorder $W$ at which the condition (17)
is satisfied decreases as the interaction strength $U$ increases. Taking
disorder value where ${\cal R}=1$ as a critical point, we obtain an estimate
for the phase diagram shown in Figure 3(b), which predicts localization for
strong interaction and strong disorder. Note, that our considerations assume a
homogeneous initial distribution of $d$-bosons, as discussed at the end of
Sec. II. Figure 3(b) suggests that at disorder $W=6.5$ system is expected to
be localized for sufficiently large $U\gtrapprox 1.5$. In what follows we
explore the properties of the localized system for $U=12$, a point that is
located deep in the localized regime according to our stability condition.

##  IV Quench dynamics as a probe of localization

The static Hartree approximation introduced above allowed to formulate a
criterion for the stability of localization. Below we proceed with probing the
dynamics of the full model in the localized regime. First, we extend the
Hartree approximation attempting to include dynamical effects in the system.
More importantly, we study the full quantum dynamics using TEBD. In this
section we fix $U=12$ and $W=6.5$ and investigate the dynamics of a quench
from the state (5) with density $\nu_{d}=1/3$.

###  IV.1 Time-dependent Hartree approximation

The static Hartree approximation relied on the assumption of quick
equilibration of the $d$-boson density. This overestimated localization in the
system, as we shall demonstrate below. Here, we adopt a different
approximation: we assume that the full many-body wave function can be
decomposed into a product of $c$\- and $d$-bosons wave functions respectively.
Such representation completely ignores entanglement between two boson species.
However, it allows performing fast and efficient simulation of dynamics on
long timescales for large systems. We note here that similar approaches to the
approximation of dynamics in MBL system have already been used [43], showing
results qualitatively similar to what we present below.

The product state structure of the wave function of the full system allows to
reduce the full many-body Schrödinger equation to the simultaneous evolution
of two non-interacting but time dependent Hamiltonians. Specifically, the time
evolution of $c$\- and $d$-bosons is governed by time-dependent Hamiltonians
$\hat{H}^{\text{dH}}_{c}(t)$ and $\hat{H}^{\text{dH}}_{d}(t)$ respectively:

$\begin{split}\hat{H}^{\text{dH}}_{d}(t)&=t_{d}\sum_{i=1}^{L-1}(\hat{d}^{\dagger}_{i+1}\hat{d}_{i}+\text{h.c.})+\sum_{i=1}^{L}V_{d,i}(t)\hat{n}_{d,i},\\\
&V_{d,i}(t)=\epsilon_{i}+U\langle\hat{n}_{c,i}(t)\rangle\end{split}$ (19)
$\begin{split}\hat{H}^{\text{dH}}_{c}(t)&=t_{c}\sum_{i=1}^{L-1}(\hat{c}^{\dagger}_{i+1}\hat{c}_{i}+\text{h.c.})+\sum_{i=1}^{L}V_{c,i}(t)\hat{n}_{c,i},\\\
&V_{c,i}(t)=U\langle\hat{n}_{d,i}(t)\rangle,\end{split}$ (20)

where the expectation value of densities of corresponding boson species reads
$\langle\hat{n}_{c/d,i}(t)\rangle=\bra{\psi_{c/d}(t)}\hat{n}_{c/d}\ket{\psi_{c/d}(t)}$,
and $\ket{\psi_{c/d}(t)}$ is the wave function obtained from the time-
evolution of the initial state (5) with the correspondent time-dependent
Hamiltonian:

$\ket{\psi_{c/d}(t)}=\hat{T}\exp\Bigr{(}-\imath\int_{0}^{t}dt^{\prime}\hat{H}^{\text{dH}}_{c/d}(t^{\prime})\Bigr{)}\ket{\psi^{0}_{c/d}}.$
(21)

Although the $d$-bosons formally are described by a many-body wave function
due to their finite density, one can exploit their non-interacting nature to
perform an efficient simulation of dynamics. In fact, as explained after Eq.
(6), both the initial state and the eigenstates can be written as product
states. After rotating the orbitals into the real space basis and back, the
time-evolution reduces to a double sum over initially occupied sites and
orbitals. The $d$-bosons can, then, be treated separately as single particles,
and their overall wave-function can be reconstructed from the single-particle
evolved states, thus reducing the complexity from exponential to linear in
$L$.

Figure 4:  Evolution of the density profile of the clean boson in the dynamic
Hartree approximation has a clear diffusive behavior. The red line corresponds
to the diffusive “lightcone”, $\langle\hat{n}_{c,x}(t)\rangle\sim\sqrt{t}$.
The data is averaged over $400$ disorder realizations, $L=90$.

First, we study the dynamics of the $c$-boson, illustrated in Figure 4, where
we plot the density profile of the clean boson
$\langle\hat{n}_{c,x}(t)\rangle$ as a function of distance from its initial
location, $x=i-L/2$, and time. We observe that in contrast to the case of
static Hartree approximation, where the $c$-boson spreads within an
exponentially localized envelope (see Appendix A), the dynamic Hartree
approximation results in a diffusive spreading of $c$-boson and its complete
delocalization over the entire system on a timescale $t\sim L^{2}$. Once the
$c$-boson spreads over the whole chain, it does not reach a steady state,
which can be attributed to the fact that this approximation oversimplifies the
true many-body character of the problem.

Figure 5:  Fluctuations of the $c$-boson density decay slightly faster than
polynomially in time. Inset shows that the saturation value of density
fluctuations reached at late times decreases monotonically with the system
size as $1/L^{b}$, with $b\approx 1.6$. This data is obtained averaging over
$400$ disorder realizations.

The fluctuations of the $c$-boson density result in a time-dependent potential
acting on disordered bosons, see Eq. (19). In order to quantify the density
fluctuations of the clean boson, we average the absolute value of the
deviation of its density from the mean value, $\delta
n_{c,i}(t)=|\langle\hat{n}_{c,i}(t)\rangle-\overline{\langle\hat{n}_{c,i}\rangle}|$,
where $\overline{\langle\hat{n}_{c,i}\rangle}$ corresponds to the average over
the interval $t\in[9.9\times 10^{3},10^{4}]$. As shown in Figure 5, the
density fluctuations around the initial position of the $c$-boson decay in
time, until they eventually reach a plateau. Note the long timescale involved
in reaching the saturation value, even for a relatively small system of $L=30$
this happens at times $t\gtrsim 10^{3}$.

The dynamics of the $d$-bosons can be intuitively understood as resulting from
an Anderson insulator globally coupled to a weak, non-Markovian, local noise.
While it is known [44, 45, 46, 47, 48] that Anderson localization is unstable
with respect to global noise, a recent work [49] has demonstrated that
coupling of an Anderson insulator to a local Markovian white noise leads to a
logarithmically slow particle transport and entanglement growth. Our dynamics
differs from that of Ref. [49] in that $d$-bosons are coupled to fluctuations
of the $c$-boson density globally throughout the entire length of the chain.
Another important difference is that the density fluctuations of the clean
boson produce a temporally and spatially correlated noise.

In Figure 6(a) we show the density profile of $d$-bosons at late times,
$T=10^{4}$, for different system sizes $L=30$, $60$, and $90$. Although the
relaxation is strongest in the middle of the chain, the memory of initial
density wave configuration survives even at late times. To understand the
dynamics of relaxation of density profile, we consider the imbalance $I(t)$
[10]. Imbalance quantifies the memory of the initial state using the
difference in occupation among initially occupied and empty sites ($N_{o}$ and
$N_{e}$ respectively),

$I(t)=\frac{N_{o}(t)-N_{e}(t)}{N_{o}(t)+N_{e}(t)}.$ (22)

For the initial period-3 density wave state considered here, the explicit form
of $N_{o/e}$ reads:

$N_{o}=\sum_{i=1}^{L/3}\hat{n}_{d,3i-2}\;,\quad
N_{e}=\frac{1}{2}\sum_{i=1}^{L/3}\Bigr{(}\hat{n}_{d,3i}+\hat{n}_{d,3i-1}\Bigr{)}.$
(23)

Figure 6(b) reveals that $I(t)$ decays without any signs of saturation even at
times $T=10^{4}$. After rescaling the time axis with a factor of $1/L^{2}$ we
observe the collapse of the data, which further supports the intuition that
although the imbalance is not a conserved quantity, its relaxation will happen
on timescale that scales as $t\propto L^{2}$ with the system size. The decay
of imbalance, after a plateau extending up to $t\approx 0.1L^{2}$ follows a
power-law scaling, $I(t)\sim\left({L^{2}}/{t}\right)^{\beta}$, with
$\beta\approx 0.4$ as shown in the inset of Fig. 8(a).

Figure 6:  (a) The density profile of $d-$bosons at late time $T=10^{4}$
retains memory of the initial state. Comparing data for different system
sizes, we observe that larger systems keep stronger memory of the initial
state. (b) Inset shows that the imbalance in the middle of the chain follows
an approximately power-law relaxation without any signs of saturation at late
times. The main panel reveals that imbalance dynamics collapses for different
system sizes after rescaling of the time axis by $1/L^{2}$. The averaging is
performed over $400$ disorder realizations.

In summary, the dynamic Hartree approximation suggests a complete
delocalization of the clean particle, that is followed by the melting of the
imbalance of disordered bosons. This does not agree with our expectations from
the criterion for stability of localization obtained in Sec. III.3, and also
disagrees with the results of the TEBD dynamics presented below. It is natural
to attribute the delocalization observed here to the nature of dynamic Hartree
approximation, that completely neglects the many-body character of the problem
and discards the correlations between the two species of bosons. The inclusion
of such correlations is expected to lead to a more efficient relaxation of
density fluctuations and may stop the diffusive spreading of $c$-boson as we
illustrate below.

###  IV.2 Stability of localization from TEBD dynamics

We now turn our attention to the fully interacting system, thus reintroducing
all the effects neglected in the Hartree approximation (both static and time-
dependent). In order to simulate dynamics of large systems we use the time-
evolving block decimation (TEBD) method [50] with time-step $dt=0.05$,
truncation $\varepsilon=10^{-9}$ and maximum bond dimension $\chi=3000$,
delegating all details of numerical implementation to Appendix B.

####  IV.2.1 Stability of localization of clean and disordered bosons

Studying the dynamics of the $c$-boson in the Hartree limit, see Appendix A,
we realized that it presents different behaviors according to the distance
from the central site. Close to the center, within a time-dependent radius
$R(t)$, the density is exponentially suppressed with a decay length
$\ell_{c}(t)$. On the other hand, farther from the middle, it shows rather a
Gaussian expansion. In the accompanying work, Ref. [36], we analyzed the
dynamical localization of the $c$-boson observed in the TEBD dynamics.
Observing features similar to the static Hartree case, there we proposed the
following ansatz for the dynamics of clean boson density,

$n_{c}(x,t)=\mathcal{N}_{c}(t)\exp\Bigr{(}-\frac{|x|}{\ell_{c}(t)\tanh\bigr{(}R(t)/|x|\bigr{)}}\Bigr{)},$
(24)

where $\mathcal{N}_{c}(x,t)$ is a normalization factor. Eq. (24) smoothly
interpolates between the two behaviors observed, as $|x|$ runs across $R(t)$.
Based on phenomenological observations, we suggested the following form for
the decay length

$\ell_{c}(t)=\xi_{c}\frac{\log\bigr{(}1+t/T_{0}\bigr{)}}{1+\log\bigr{(}1+t/T_{0}\bigr{)}},$
(25)

which implies saturation of the decay length to a constant localization length
$\xi_{c}$ on timescales $t\gg T_{0}$.

Figure 7:  The spreading of the $c$-boson density in the $L=60$ system is well
described by the ansatz of Eq. (24) (blue curve), while it deviates from the
diffusive behavior (red curve).

Figure 7 shows the evolution of the density profile of the clean boson with
time. While at early times it approximately agrees with a diffusive profile
extracted from time-dependent Hartree approximation (red line and Fig. 4), at
times of order one, the density spreading slows down. To describe the envelope
of the density profile, we use the fact that when $R(t)$ saturates to $L/2$,
the ansatz (24) reduces to $n_{c}(x,t)\approx e^{-|x|/\ell_{c}(t)}$. By
imposing the condition $n_{c}(x,t)=\text{const}$, we obtain the scaling
$|x|\approx\ell_{c}(t)$ shown in Fig. 7 by the blue line, that provides better
description for the envelope of the density profile.

Figure 8:  (a) TEBD dynamics (blue line) reveals a broad plateau in imbalance
at late times, whereas time-dependent Hartree approximation (red line)
predicts a power-law decay. (b) The imbalance of central six sites in TEBD
dynamics reveals a weak power-law decay. The stronger deviation between TEBD
and approximate time-dependent Hartree data highlights the importance of
entanglement for the accurate description of the dynamics. Data is shown for
$L=60$ and averaged over $400$ and $50$ disorder realizations for time-
dependent Hartree and TEBD respectively.

Next, we consider the dynamics of imbalance of $d$-bosons, defined in Eq.
(22), that is shown in Fig. 8(a). While TEBD and time-dependent Hartree
approximation approximately agree at times up to $t\approx 10$, at later times
the imbalance obtained from TEBD remains approximately constant, while the
time-dependent Hartree approximation predicts a power-law decay of the
imbalance. This is further highlighted by the exponent $\beta$, shown in the
inset, that is compatible within error bars with value $\beta=0$ for TEBD
dynamics, which corresponds to non-decaying imbalance.

In order to quantify the stronger relaxation close to the center of the chain,
triggered by the presence of the clean boson, we also study the local
imbalance in the middle of the chain, $I_{\text{mid}}(t)$ defined according to
Eq. (22) but restraining the sum in Eq. (23) to the six central sites. The
quantity $I_{\text{mid}}(t)$ shows a stronger decay in Fig. 8(b) as compared
to its global counterpart in Fig. 8(a). However, again the TEBD dynamics
reveals a much weaker effect of the $c$-boson compared to the Hartree
approximation. Fitting the decay to the power-law form, the resulting exponent
remains small and it does not change with system size (inset). Hence, we
expect the effect on the central region to be system size independent, thus
leaving the boundaries unaffected as $L\to\infty$.

####  IV.2.2 Entanglement dynamics

After demonstrating the stability of localization in TEBD dynamics, we present
additional details for the picture of entanglement spreading provided in Ref.
[36]. Thanks to $U(1)$ conservation, we can write the reduced density matrix
of the first $i$ sites from the left, $\rho_{i}$, in block-diagonal form,
$\rho_{i}=\sum_{n}p_{i}^{(n)}\rho_{i}^{(n)}$, where $p_{i}^{(n)}$ corresponds
to the weight of the $n$-particles sector and $\mathrm{tr}\rho_{i}^{(n)}=1$.
It is then convenient to separate the entanglement entropy,
$S(i)=-\mathrm{tr}\rho_{i}\log\rho_{i}$, in two contributions: the
configuration entanglement $S_{C}(i)$ and the particle entanglement $S_{n}(i)$
[51]

$\displaystyle S(i)=$ $\displaystyle S_{n}(i)+S_{C}(i)$ (26) $\displaystyle
S_{C}(i)=-$
$\displaystyle\sum_{n}p_{i}^{(n)}\mathrm{tr}\rho_{i}^{(n)}\log\rho_{i}^{(n)}$
(27) $\displaystyle S_{n}(i)=$ $\displaystyle-\sum_{n}p_{i}^{(n)}\log
p_{i}^{(n)}.$ (28)

Figure 9:  (a) The configuration entropy shows a steady logarithmic growth,
whose start is delayed as the cut gets farther from the center of the chain,
due to the weaker effect of the interaction at the boundaries. (b) The real
space profile of the configuration entropy is highly non-uniform, with entropy
being maximal in the center of the chain, where the $c$-boson is localized.
Far from the center, $S_{C}(i)$ is close to zero, indicating the absence of
interactions in the boundaries up to late times. (c) The growth of number
entropy is similarly delayed with increasing distance from the center of the
chain. It shows slow growth, compatible with the weak relaxation of $d$-bosons
observed in the imbalance dynamics. (d) The number entropy also has a non-
uniform profile, with its values far from the center being consistent with the
particle entropy of the Anderson insulator. This data refer to a system of
$L=60$ sites, averaged over $50$ disorder realizations.  Figure 10:  The
configuration entropy $S_{C}(i,t)$ at $U=12$ and $L=60$, characterizing the
interacting nature of the system, presents an extremely inhomogeneous growth.
In agreement with the picture of propagation of MBL, its growth is constrained
within the “MBL light-front” $x_{\text{MBL}}(t)$ shown by the red line. When
$x\gg x_{\text{MBL}}$, the configuration entropy $S_{C}$ is close to zero,
suggesting that system remains effectively non-interacting until late times.

The configuration entropy $S_{C}(i,t)$ arises from the superposition of
different product states in $\ket{\psi(t)}$ and measures the correlation
between particle configurations in the two subsystems. The growths of
configuration entropy can be attributed to the interacting nature of the
system. Thus configuration entropy does not show interesting dynamics in
Anderson insulator, while displaying a logarithmic growth in a many-body
localized phase [51] providing the main contribution to the entanglement
growth [6]. Particle number entanglement entropy $S_{n}(i)$ arises from the
occupation of different sub-sectors of the density matrix $\rho_{i}$ and hence
accounts for the particle transport in the system. Recent studies [52, 53]
suggested an extremely slow, albeit finite, growth $S_{n}(i,t)\sim\log\log t$,
which they attributed to slow particle transport. Subsequent work by Luitz and
Lev [54] however suggested absence of such transport and saturating particle
number entanglement.

The dynamics of the two contributions to the entanglement in our model is
presented in Figure 9, where the first row shows the configuration
entanglement growth for different cuts in the chain (a) and its profile at
different times (b). In the second row, we show the dynamics of particle
number entanglement. Similarly to the configuration entropy, $S_{n}(i,t)$ has
a non-homogeneous growth across the system shown in Fig. 9(c)-(d).

Separating entanglement $S(i)$ into configurational and particle part, it is
clear that the logarithmic growth is due to the configuration entropy,
$S_{C}(i)\sim\xi_{S}\log(t)$, as clearly shown in Figure 9(a). The
configuration entanglement growth has a non-homogeneous behavior in the
system, explained by the propagation of MBL phenomenology discussed in Ref.
[36]. Introducing an effective interaction $U_{\text{eff}}(x,t)=Un_{c}(x,t)$,
the dynamics of $S_{C}(i)$ can be described by the ansatz
$S_{C}(i,t)\approx\xi_{S}\log(1+U_{\text{eff}}(i,t)t)$. This explains the
delay of onset of growth of configurational entropy in Fig. 9(a) away from the
center of the chain. In addition, this phenomenology predicts an linearly
decreasing entanglement profile away from the center of the chain. Indeed,
using an exponential profile for the density of the clean boson at late times
we obtain the following effective interaction strength, $U_{\text{eff}}\approx
Ue^{-|i-L/2|/\ell_{c}(t)}$. This gives the following asymptotic behavior of
configurational entanglement,

$S_{C}(i,t\gg 1)\sim\text{const}-\frac{\xi_{S}}{\ell_{c}(t)}|i-L/2|,$ (29)

that decays linearly away from the center of the chain with slope given by
$1/\xi_{S}$. This prediction is confirmed by the numerical results shown in
Fig. 9(b), yielding the value of $\xi_{S}\approx 0.31$.

According to the effective interaction picture describing the propagation of
MBL, the interacting nature spreads through the system producing a many-body
localization lightcone. At late times, when
$n_{c}(i,t)\approx\mathcal{N}_{c}(t)e^{-|i-L/2|/\ell_{c}(t)}$, this can be
captured analytically by solving $S_{C}(x,t)=\text{const}$ for $x$ and defines
the MBL “light-front”,

$x_{\text{MBL}}(t)=\ell_{c}(t)\log\bigr{(}\mathcal{N}_{c}(t)Ut\bigr{)}.$ (30)

As shown in Figure 10, the prediction of the MBL lightcone (red curve)
provides an accurate description of the actual behavior of the spread of
configuration entropy, thus testing the genuine propagation of the many-body
nature through the chain.

##  V Probing eigenstates with DMRG-X

Although TEBD enabled simulation of dynamics for large systems, the maximal
time remained limited by entanglement growth. Below we use the DMRG-X
algorithm [38, 39, 40] to probe eigenstates of large systems, providing an
effective insight into infinite time behavior.

Figure 11:  The density profile of the $c$-boson averaged over different
eigenstates as a function of the distance from $i_{\text{max}}$ shows an
exponential decay for both the system sizes considered. Through a fit to the
density in different regions we obtain an upper and lower bound for the
localization length, with $\xi^{\text{low}}_{c}\approx 2.1$ (red dashed line)
corresponding to the fit excluding the central site and
$\xi^{\text{up}}_{c}\approx 3$ (black dashed line) obtained from the fit that
exclude 12 central sites.

###  V.1 $c$\- and $d$-boson localization

We extract highly excited eigenstates of the Hamiltonian (4) with $L=30$ and
$60$ sites, $U=12$ and 1/3 filling of $d$-bosons. We obtain in total $250$
(500) eigenstates in the middle of the spectrum by targeting $25$ different
energies for each of the $10$ ($20$) disorder realizations considered for
$L=30$ (60) chain respectively. 111We note, that the first 10 disorder
realizations for $L=60$ system are chosen in such a way that the values of the
random energies of the $30$ central sites agree with those for $L=30$ system.
The states are obtained as a result of $100$ sweeps where the targeted
eigenstate is variationally approximated by an MPS of maximum bond dimension
$\chi=500$ and $250$ for systems of size $30$ and $60$, respectively. The
algorithm is initialized with different copies of the initial state (5), where
the $c$-boson initial position is chosen among empty sites in the central
region. We refer the reader to Appendix C for details on the implementation of
the algorithm as well as its performance metrics.

For each of the eigenstates we extract the peak position of the $c$-boson
density, $i_{\text{max}}$. The average $c$-boson density plotted as a function
of the distance from its peak presents an exponential profile, as shown in
Fig. 11. We notice that the density profile around the peak is enhanced with
respect to the tails, which can be potentially attributed to the large value
of $U$ producing stable doublons. The exponentially decaying density confirms
the localization of the clean particle and allows us to extract the upper and
lower bound on the localization length, $\xi^{\text{low}}_{c}\approx 2.1$ and
$\xi^{\text{up}}_{c}\approx 3$ resulting from fitting its profile close to the
center (excluding the central site) and at the edges of the chain. The two
bounds quantitatively agree with the lengthscale
$\ell_{c}(t\to\infty)=\xi_{c}=2.5$ _extrapolated_ from TEBD dynamics in Ref.
[36], suggesting that the localization of the $c$-boson observed in dynamics
is not a transient effect, but it is a property of the system, provided the
correct sector of the Hilbert space is explored.

Figure 12:  Strong deviation of the $d$-bosons density away from the thermal
value $\nu_{d}=1/3$ consistent between different system sizes indicates the
breakdown of thermalization in eigenstates.

Probing potential localization of $d$-bosons in eigenstates is more
complicated: the finite particle density does not allow to define a
localization center. Therefore, we consider the average deviations of the
$d$-bosons density from the thermal value given by their average density
$\nu_{d}=1/3$ as a function of the distance from $i_{\text{max}}$. In Figure
12, we show the averaged absolute value of the deviation of the density
expectation value from $\nu_{d}$, which is limited to the range of values
$0<\langle|\hat{n}_{d,i}-\nu|\rangle<1-\nu_{d}$. While in a thermalizing
system we expect this quantity to be small and decrease exponentially with the
system size, the eigenstates of our problem have on average large expectation
value of $\langle|\hat{n}_{d,i}-\nu|\rangle$. In addition, increasing the
system size does not lead to a decrease of the average distance from the
thermal value, thus suggesting that the $d$-bosons density strongly fluctuates
around $\nu_{d}$ in eigenstates due to their localization. Interestingly, we
notice a slight peak of this quantity at $i=i_{\text{max}}$, which is in
agreement with the enhancement of the $c$-boson density observed in Fig. 11
and may be attributed to the effect of doublons. In the region around the
center, however, we observe a slight weakening of localization, corresponding
to smaller values of $\langle|\hat{n}_{d,i}-\nu|\rangle$. This can be
attributed to the fact that in the vicinity of the peak of the $c$-boson, the
interactions are effectively stronger, leading to an increased relaxation of
the $d$-bosons.

Figure 13:  The spectrum of $\rho_{d}$ shows a step behavior, dropping quickly
to small values at $i/L=\nu_{d}=1/3$ (dashed line). This corresponds to the
presence of $N_{d}$ occupied single-particle orbitals, thus confirming
localization. We note that $w_{d,i}>1$ is due to the bosonic nature of these
particles. The inset shows the log-averaged spectrum of $\rho_{c}$ on a
logarithmic scale, revealing that there is a single significantly occupied
orbital, while all the others have an exponentially decaying weight.  Figure
14:  (a) Average entanglement as a function of the cut position results in a
featureless profile. The agreement in entanglement between different system
sizes where they overlap is due to particular disorder choice. (b) Averaging
the entanglement profile as a function of the distance from the center of
localization of $c$-boson, $i_{\text{max}}$, shows enhancement in entanglement
caused by the presence of the clean boson. The dashed line shows the
comparison to an exponential fit. The collapse of entanglement profile between
different system sizes suggests area-law entanglement scaling.

Further evidence of localization for both types of particles is found in the
spectrum of the single-particle density matrices
$\rho_{c}=\langle\hat{c}^{\dagger}_{i}\hat{c}_{j}\rangle$
($\rho_{d}=\langle\hat{d}^{\dagger}_{i}\hat{d}_{j}\rangle$) [55, 56]. In the
particular case treated here, the von Neumann entropy of $\rho_{c}$
corresponds to the intra-species entanglement
$S_{cd}=-\operatorname{tr}\rho_{c}\log\rho_{c}$. We study its distribution and
average value among eigenstates, obtaining average entropy of $S_{cd}\approx
0.85$ for both system sizes. Furthermore, the eigenvalues of $\rho_{c/d}$,
$w_{c/d,i}$, sorted by decreasing value, can be interpreted as occupation
numbers of single-particle orbitals [55, 56]. In the case of an MBL system, it
is expected that particles sit on almost localized sites, hence only the first
$N_{c/d}$ orbitals should be significantly occupied. In Figure 13 we show the
average, $\langle\cdot\rangle$, and log-average, $\overline{x}=e^{\langle\ln
x\rangle}$, value of the ordered eigenvalues for $\rho_{d}$ and $\rho_{c}$
respectively. The scaling of $\langle w_{d,i}\rangle$ shown in the inset
presents a step behavior, consistent with the localization of the $d$-bosons.
Similarly, the typical value $\overline{w}_{c,i}$ presents a single large
eigenvalue, while the rest decay exponentially.

###  V.2 Entanglement structure of eigenstates

In Figure 14(a) we show the bipartite entanglement profile averaged as a
function of the position of the cut $i$. While the entanglement profile of
eigenstates suggests an area-law entanglement scaling characteristic of many-
body localization [7, 57], this representation of data does not show any
effect from the presence of the clean boson. To observe the influence of the
$c$-boson, we average the entanglement profile defined with respect to the
distance from the localization site of the $c$-boson, $i-i_{\text{max}}$.
Figure 14(b) reveals a peak in $S(i-i_{\text{max}})$ around zero, indicating
that the clean boson is responsible for additional entanglement in
eigenstates. Away from $i_{\text{max}}$ entanglement saturates to a constant
value, revealing a clear area-law scaling. Fitting the decay of entanglement
away from the center with an exponential of the type
$S_{0}+ce^{-|i-i_{\text{max}}|/\zeta_{S}}$ allows to extract the value of
$\zeta_{S}\approx 4$, suggesting that the single $c$-boson is capable of
generating non-trivial entanglement patterns on a scale larger than its
localization length.

The increase in entanglement of eigenstates that decays exponentially with the
distance away from the location of $c$-boson provides further support to the
localized nature of eigenstates. Moreover, this provides a complementary view
on the picture of the inhomogeneous effective interaction triggered by the
presence of $c$-boson and resulting in the dynamical “propagation of MBL”
presented in [36].

##  VI Discussion

We investigated the effect of the coupling of a small local bath represented
by a single clean boson to a one dimensional Anderson insulator with finite
particle density. First, we used a static Hartree approximation that
overestimates localization to obtain the parameters range most favorable for
localization. In addition, we analyzed the effect of two-particle resonances
triggered by the interaction, estimating the region of parameters where
localization is perturbatively stable with respect to such resonances.

Focusing on the parameter values where the perturbative criterion predicts
localization, we studied the dynamics in the coupled system. We demonstrated
that the time-dependent Hartree approximation, which neglects the entanglement
between the localized particles and the degree of freedom representing the
bath, predicts a slow delocalization of the system consistent with diffusion.
In contrast, the TEBD numerical simulation, which takes into account all
quantum correlations among the two bosonic species, shows evidence of
localization within the experimentally relevant timescales achieved. In
particular, the imbalance of disordered particles saturates, suggesting that
memory of the initial density-wave configuration is retained even at long
times. In addition simulations reveal logarithmic growth of configuration
entanglement, in agreement with the MBL phenomenology, and extremely slow
growth of particle number entanglement, ascribed to the weak relaxation of the
disordered particles due to the interaction with the small bath. In agreement
with Ref. [36], we observed an enhanced effect of the bath on the $d$-bosons
in the center of the system, while far from it the difference with the
Anderson case is negligible.

Finally, in order to probe the fate of localization at infinite times, we used
DMRG-X method to extract highly excited eigenstates of the Hamiltonian. We
observe that the clean boson remains localized in eigenstates, confirming the
intuition gained from the study of dynamics. We further notice area-law
scaling of the entanglement entropy of eigenstates, with an enhancement in the
vicinity of the localization site of the $c$-boson, $i_{\text{max}}$. These
results provide evidence for the effective stability of localization even at
infinite times and further support for the picture proposed in the
accompanying work [36].

In summary, our work provides evidence for MBL proximity effect and a
phenomenological picture of dynamics [36] for large systems that are beyond
the reach of numerical exact diagonalization. Our predictions are readily
verifiable in state of the art experimental setups [33, 34] that are capable
of probing the particle dynamics with a single-site resolution.

However, there remain many open questions. Describing a more “powerful”
thermal bath requires a finite $c$-bosons density. The persistence of the
phenomenology proposed here in such a scenario is non-trivial and its study
can lead to a better understanding of the MBL proximity effect. In the same
direction, a careful investigation of the role of the interaction strength is
in order for an accurate description of MBL-bath systems. In particular, in
the opposite limit of weak interaction the $c$-boson is expected to delocalize
quickly yielding a very weak and uniform effective coupling. This regime can
provide a different point of view on the MBL-bath model. A further interesting
direction would be to address the stability of our predictions to inter-
species interactions, i.e. replacing the Anderson insulator Hamiltonian
$\hat{H}_{d}$ with an MBL Hamiltonian. Although we expect the qualitative
picture discussed above to hold, a different entanglement dynamics is supposed
to arise in this case. Finally, replacing the Hamiltonian evolution with a
periodic Floquet drive may allow the investigation of the long-time behavior
of the model, since it would avoid the large number of singular value
decompositions required for reaching long times with TEBD algorithm. In this
framework, the saturation of the measured quantities could be observed, thus
providing insight into the infinite-time fate of the system.

###### Acknowledgements.

We thank M. Ljubotina for insightful discussions. P. B., A. M. and M. S.
acknowledge support by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (Grant Agreement No.
850899). D. A. was supported by the Swiss National Science Foundation and by
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (Grant Agreement No. 864597). The development
of parallel TEBD code was supported by S. Elefante from the Scientific
Computing (SciComp) that is part of Scientific Service Units (SSU) of IST
Austria. Some of the computations were performed on the Baobab cluster of the
University of Geneva.

##  Appendix A Clean boson dynamics in static Hartree disorder

Figure 15:  Hartree induced localization length compared with the analytic
result of Eq.(9), dashed lines, reveals good agreement at weak disorder, while
at larger $W$ the deviation of $P(\tilde{\epsilon}_{i})$ from Gaussian leads
to a worse agreement. The localization length $\xi_{c}$ reflects the behavior
of the effective disorder, leading to larger $\xi_{c}$ for systems far from
half filling. The weaker nature of the effective disorder results in
$\xi_{c}\gg\xi_{d}$ for all the parameters studied, thus highlighting the much
weaker localization of the $c$-boson.

In the main part of the current paper, we discussed the static Hartree
approximation, yielding an effective disorder acting on the $c$-boson and
consequently its localization. As the effective disorder changes with the
filling fraction of the $d$-bosons, we expect also the localization length to
change accordingly. This should lead to a weaker localization whenever the
$d$-bosons density, $\nu_{d}$, deviates from half filling. This is confirmed
by the numerical results shown in Figure 15, where we observe stronger
localization for $\nu_{d}=1/2$. Additionally, we compare the data with the
prediction of Eq. (9) in the main text (dashed lines). Although qualitatively
similar, the quantitative agreement fails as the disorder strength $W$
increases. We attribute this discrepancy to the deviation of the distribution
of effective disorder from Gaussian discussed in the main text.

Figure 16:  The density profiles (a) of the single particle in the static
Hartree disorder at $U=2.5$ show an exponential decay within a range $R(t)$,
outside of which they rather have a Gaussian behavior. The exponential decay
is characterized by a decay length $\ell_{c}(t)$ (b), which after an initial
transient is constant in time and saturates to the value obtained in Fig. 15.
The extent of the exponentially suppressed region, $R(t)$, instead grows in a
power-law fashion until it reaches the boundaries of the system (c).

Next, we consider the $c$-boson dynamics in the static Hartree approximation
fixing the disorder strength to $W=6.5$, close to the minimum of its
localization length. The dynamics in this setup can be obtained numerically
for large systems without resorting to large computational resources, as the
single-particle Hilbert space dimension only scales as the system size $L$.
Therefore, we study the density dynamics of a $c$-boson initialized in the
center of a $L=120$ chain evolving the initial state up to $t\approx 10^{4}$.
The density profiles obtained in this way in Figure 16(a) reveal the presence
of two distinct behaviors, depending on the distance from the center. Close to
the initial position of the particle, the density in Fig. 16(a) decays
exponentially, consistently with the expectation for a localized particle.
However, further away from the initial location of the particle, the density
dependence is consistent with a Gaussian profile. While the Gaussian profile
is expected from the diffusive spreading of the particle that could take place
at early times for weak disorder (weak localization), the coherent
backscattering from static disorder leads to the localization of the particle.

The two distinct regimes in the density profile can be described by
introducing two separate lengthscales: the distance $R(t)$ corresponds to
crossover from exponential to diffusive behavior, and $\ell_{c}(t)$ describes
the slope of exponential decay. Using these lengthscales, the dynamics of the
density can be captured by the following ansatz,

$n(x,t)\propto\exp\Bigr{(}-\frac{|x|}{\ell_{c}(t)\tanh\bigr{(}R(t)/|x|\bigr{)}}\Bigr{)},$
(31)

that smoothly interpolates between exponential and Gaussian decays for
$x<R(t)$ and $x>R(t)$ respectively.

The ansatz (31) can then be used to fit the numerical data and to obtain the
behavior of the two parameters $\ell_{c}(t)$ and $R(t)$, shown in Figure
16(b)-(c). We observe that for the static Hartree approximation, the decay
length after a transient quickly saturates to the corresponding localization
length $\xi_{c}$ shown in Figure 15 by a dashed line. The extent of the
exponential decay region around the center, $R(t)$, instead increases as a
power-law, before its saturation to half system size, meaning that the
localized region extends over the whole chain. In the companion paper [36] and
in the main text we use the same ansatz to study the dynamics in the fully
interacting system, obtaining similar results, although interactions result in
much more intricate dynamics.

##  Appendix B Details and benchmarks of TEBD implementation

Figure 17:  Cartoon representation of the TEBD algorithm. The time evolution
can be naturally parallelized, by sharing each of the operations contoured in
red to a different core.

We use the TEBD algorithm [50] to perform the time evolution of the wave
function of the large system represented as an MPS. The algorithm splits the
unitary evolution into time steps (Trotterization), $dt$, and further divides
them into even and odd gates, $\hat{G}_{i,i+1}=\exp{\bigr{(}-\imath
dt\hat{h}_{i,i+1}\bigr{)}}$, where $\hat{h}_{i,i+1}$ is the Hamiltonian
density on sites $i$ and $i+1$. Since even(odd) gates do not overlap with one
another, all even (odd) gates can be applied independently to the state. This
allows a natural parallel implementation of the algorithm, where all even
(odd) gates are applied simultaneously.

The TEBD algorithm has two main sources of error: the finite value of time-
step used in Trotterization and the truncation error in the singular-value
decomposition (SVD). The time step error naïvely scales as $O(dt^{2})$, but it
can be reduced to $O(dt^{4})$ with an appropriate higher order Suzuki-Trotter
decomposition. By choosing sufficiently small time steps, then, the accuracy
of the simulation remains reliable for long timescales. The accuracy of the
simulation can be controlled by the degree of energy conservation violation,
i.e. by tracking $\Delta
E(t)=|\bra{\psi_{0}}\hat{H}\ket{\psi_{0}}-\bra{\psi(t)}\hat{H}\ket{\psi(t)}$.
In our simulations, we fixed $dt=0.05$. In spite of a slower execution time,
such a small time-step ensures very accurate results, as evidenced by $\Delta
E(T)<10^{-4}$ at the final simulation time $T=200$.

The second source of error, i.e. the singular values truncation $\varepsilon$,
arises from the singular-value decomposition performed after each gates
application in order to restore the MPS form. In order to avoid the
exponential growth of the MPS bond dimension $\chi$, one needs to select the
most relevant singular values by neglecting all those which are smaller than a
certain threshold $\varepsilon$. The truncation then corresponds to neglecting
some weight of the wave-function, and the error can be estimated by the sum of
the discarded singular values $\sum_{\lambda_{i}<\varepsilon}\lambda_{i}^{2}$.
A further limit on the number of singular values is imposed by the maximal
bond dimension, $\chi$. If during the time evolution the number of singular
values smaller than $\varepsilon$ becomes larger than $\chi$, i.e. the bond
dimension is saturated, then the control on the accuracy of the results is
lost, as there is no guarantee that large singular values will be discarded in
the future. In our simulations, we fix $\epsilon=10^{-9}$ and $\chi=3000$.
This set of parameters guarantees the reliability of the results as the
maximum bond dimension is saturated only in the last 10 time-steps (in hopping
parameter units) of less than $1/10$ of the central bonds. Furthermore, we can
ensure accuracy by comparison to exact diagonalization results on smaller
systems, which for the chosen simulation parameters show discrepancies in
local observables and entanglement typically of the order $O(10^{-4})$.

The choice of these parameters, however, makes the simulations extremely
demanding. However, the structure of the algorithm is naturally suited for
parallelization, as shown in Figure 17, and the use of a parallel
implementation of the code allowed us to study the dynamics even in the high
entanglement regime, where most of the tensors have a large bond dimension of
order $1000$.

##  Appendix C Details and benchmarks of DMRG-X implementation

The DMRG-X algorithm [38, 39, 40] used to obtain the highly excited
eigenstates of the Hamiltonian (4) relies on the shift-invert method applied
to MPS states. The variant of the algorithm used in this work is described in
detail in the supplementary material of [38]. In this section, we will give a
brief overview of the algorithm and then explore the dependence of the quality
of the results with respect to various tuning parameters such as bond
dimension, system size, etc.

The DMRG-X algorithm used can be understood as a simple modification of the
standard two-site DMRG [58] algorithm which is extensively used to target
ground states of one-dimensional systems. We describe the two differences
compared to the ground state algorithm.

First, the ground state DMRG uses the matrix product operator (MPO) form of
the Hamiltonian, $H$, to perform the variational optimization. The DMRG-X
algorithm, in contrary, employs $\mathcal{H}=(H-E_{0})^{2}$, which is still
efficiently represented by a MPO with a slightly larger bond dimension. The
energy $E_{0}$ is the target energy, i.e. the energy of the eigenstates we are
aiming to calculate. We square the operator in order to make it semi-positive
definite. This is done in order to use iterative inversion algorithms (see
below) which require such property and can be avoided if a different iterative
algorithm is chosen. We have numerically observed that different iterative
algorithms are less stable and converge slower for most cases.

Figure 18:  (a) Increasing the bond dimension $\chi$ leads to a better
distribution of the variance of the energy of the states $\sigma^{2}$,
shifting it to lower values. Similarly, the increase of the number of
conjugate gradient iterations at fixed $\chi$ improves the dataset quality.
(b) The distribution of the distance of the final energy $E$ from the target
energy $E_{0}$. At fixed bond dimension $\chi=250$ doubling the number of
conjugate gradient iterations produces a dramatic improvement of the final
energy distribution. Further increasing the bond dimension to $\chi=500$ then
gives an additional refinement of the distribution. (c) The entanglement
distribution $P(S)$ for the bonds closest to $i_{\text{max}}$ becomes narrower
as the number of conjugate gradient iterations is increased. However, as the
bond dimension is increased at fixed $CG$ the MPS states can host larger
entanglement, as is confirmed by the broader $P(S_{max})$ at $\chi=500$.

Second, in the DMRG-X algorithm, we are not aiming to find the eigenstate of
$\mathcal{H}$ with the largest magnitude eigenvalue. Instead, the goal is to
find the eigenstate corresponding to the lowest magnitude eigenvalue, i.e. the
one which corresponds to the eigenstate closest to energy $E_{0}$. This is
taken into account in the local optimization step, present in DMRG-type
algorithms, i.e. when two adjacent local tensors are optimized. In this step,
one generates an effective operator by tracing the left and right
environments, $\mathcal{H}_{\text{eff}}$ [58]. Calculating the lowest
magnitude eigenvalue of $\mathcal{H}_{\text{eff}}$ is equivalent to
calculating the largest eigenvalue of $\mathcal{H}_{\text{eff}}^{-1}$. We use
a power-method to converge to the largest magnitude eigenvalue geometrically.
Instead of applying $\mathcal{H}_{\text{eff}}^{-1}$ at each iteration, we
solve the linear system
$\mathcal{H}_{\text{eff}}\textbf{x}_{i+1}=\textbf{x}_{i}$ using a conjugate
gradient algorithm. Due to the sweeping nature of the DMRG-type algorithms we
have observed that a single iteration of the local power method leads to a
faster convergence of the global wavefunction (for fixed runtime), indicating
that increasing the number of sweeps is a more important factor. The iterative
algorithm scales as $O(CG\chi^{3})$, where $CG$ is the number of conjugate
gradient iterations, instead of the $O(\chi^{6})$ scaling of the exact
diagonalization of $\mathcal{H}_{\text{eff}}$. This scaling, together with the
explicit conservation of both particle numbers, allow us to reach bond
dimensions up to $\chi=500$.

For the simulations presented in this work we used the following parameters:
In the global level we restricted to 100 DMRG sweeps. However, we do not
always keep the final state but the state with the least energy variance as
measured at the end of every sweep for the last 50 sweeps. This is because the
algorithm tends to switch between different eigenstates as it converges
towards the target energy and during the switching period it is not always
accurately converged. This behaviour is not present in the standard DMRG
algorithm and is attributed to the exponentially large density of states
around the target energy. For the solution of each local optimization problem,
we use a single power method iteration which is performed using the conjugate
gradient algorithm. For the singular value decomposition which is used to
extract the updated local tensors we have used a very low probability
truncation cutoff $\epsilon=10^{-28}$, in order to avoid any statistical bias
towards eigenstates of low entanglement. We found that the most important
tuning parameters are the bond dimension $\chi$ and the number of conjugate
gradient iterations $CG$. Below we show various benchmarks for different
parameters and system sizes.

We compare the quality of the resulting eigenstates for $\chi=250,500$ and
$CG=600,1200$. A first check is obtained by comparison of the energy variance

$\sigma^{2}=\langle\hat{H}^{2}\rangle-\langle\hat{H}\rangle^{2}$ (32)

of the resulting states. As true eigenstates have exactly $0$ energy variance,
a small $\sigma^{2}$ is an indicator of good convergence. In Figure 18(a), we
compare the variance distribution for the different range of parameters,
observing, as expected, that increasing both the number of iterations and the
bond dimension leads to a better quality result.

The low-valued distribution of the variance suggests that the final states
obtained with the DMRG-X algorithm are actually extremely close to exact
eigenstates of $\hat{H}$. As a next step, one must check whether the obtained
states are indeed highly excited. We thus study the distribution of the energy
difference among the obtained state and the target energy $E_{0}$. As shown in
Figure 18(b), increasing $\chi$ and $CG$ yields a narrower distribution close
to $0$.

Figure 19:  The distribution of the peak of the $c$-boson is almost
homogeneous through the whole system in the $L=30$ case, but as the system
size is doubled the typical position of $i_{\text{max}}$ is more restrained to
the center of the chain. The broad distribution of $i_{\text{max}}$ is
expected, as the local optimization performed in the algorithm can indeed
update the position of the $c$-boson from its initial value and “move” it
through the system.

Since the amount of entanglement entropy an MPS state can host is strictly
related to its bond dimension, it is important to observe how the entanglement
varies as we increase $\chi$. We thus study the distribution of the
entanglement entropy on the $10$ bonds closer to the peak of the $c$-boson,
$i_{\text{max}}$, where entanglement is supposed to be largest. The results
are shown in Figure 18(c) and show that at fixed bond dimension, a larger
number of iterations leads to a narrower distribution, in agreement with the
small entanglement of exact eigenvalues of MBL Hamiltonians. On the other
hand, fixing $CG=1200$ and increasing the bond dimension gives access to
states with larger entanglement, as verified by the broader distribution of
$S$ at $\chi=500$.

Finally, we study the distribution of the position of the peak of the
$c$-boson. As initial state in the DMRG-X algorithm we always used a version
of $\ket{\psi_{0}}$, Eq. (5), with the $c$-boson shifted of a few sites from
the center, provided the site is not already occupied by a $d$-boson to avoid
stable doublons. However, the local optimization performed by the algorithm
can easily “move” the $c$-boson, as can be seen by the almost even
distribution of $i_{\text{max}}$ shown in Figure 19, for $L=30$. This behavior
is consistent in all the datasets studied and does not seem to have a
particular relation with the quality of the eigenstates. As the system size is
doubled, the distribution of the peak of the $c$-boson becomes more centered
around $L/2$, probably due to the larger density of state which increases the
possibility of finding an eigenstate close to the target energy for smaller
displacement of the $c$-boson.

## References

  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Physical Review 109, 1492 (1958).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-$t$ transport, Phys. Rev. Lett. 95, 206603 (2005).
  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal-insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 321, 1126 (2006).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013a).
  * Serbyn  _et al._ [2013b] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013b).
  * Huse  _et al._ [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phenomenology of fully many-body-localized systems, Phys. Rev. B 90, 174202 (2014).
  * Imbrie [2016] J. Z. Imbrie, On Many-Body Localization for Quantum Spin Chains, Journal of Statistical Physics 163, 998 (2016).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Observation of many-body localization of interacting fermions in a quasirandom optical lattice, Science 349, 842 (2015).
  * Lüschen  _et al._ [2017a] H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017a).
  * Lukin  _et al._ [2019a] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Probing entanglement in a many-body–localized system, Science 364, 256 (2019a).
  * Rispoli  _et al._ [2019] M. Rispoli, A. Lukin, R. Schittko, K. Sooshin, M. E. Tai, J. Léonard, and M. Greiner, Quantum critical behaviour at the many-body localization transition, Nature 573, 385 (2019).
  * Choi  _et al._ [2016] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Exploring the many-body localization transition in two dimensions, Science 352, 1547 (2016).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Probing slow relaxation and many-body localization in two-dimensional quasiperiodic systems, Phys. Rev. X 7, 041047 (2017).
  * Guo  _et al._ [2019] Q. Guo, C. Cheng, Z.-H. Sun, Z. Song, H. Li, Z. Wang, W. Ren, H. Dong, D. Zheng, Y.-R. Zhang, R. Mondaini, H. Fan, and H. Wang, Observation of energy resolved many-body localization, arXiv e-prints , arXiv:1912.02818 (2019), arXiv:1912.02818 [quant-ph] .
  * De Roeck  _et al._ [2016] W. De Roeck, F. Huveneers, M. Müller, and M. Schiulaz, Absence of many-body mobility edges, Phys. Rev. B 93, 014203 (2016).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field heisenberg chain, Phys. Rev. B 91, 081103 (2015).
  * Brighi  _et al._ [2020] P. Brighi, D. A. Abanin, and M. Serbyn, Stability of mobility edges in disordered interacting systems, Phys. Rev. B 102, 060202 (2020).
  * Yao and Zakrzewski [2020] R. Yao and J. Zakrzewski, Many-body localization in the bose-hubbard model: Evidence for mobility edge, Phys. Rev. B 102, 014310 (2020).
  * Pomata  _et al._ [2020] N. Pomata, S. Ganeshan, and T.-C. Wei, In search of a many-body mobility edge with matrix product states in a Generalized Aubry-André model with interactions, arXiv e-prints , arXiv:2012.09853 (2020), arXiv:2012.09853 [cond-mat.dis-nn] .
  * De Roeck and Huveneers [2017] W. De Roeck and F. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Potirniche  _et al._ [2019] I.-D. Potirniche, S. Banerjee, and E. Altman, Exploration of the stability of many-body localization in $d>1$, Phys. Rev. B 99, 205149 (2019).
  * Wahl  _et al._ [2019] T. B. Wahl, A. Pal, and S. H. Simon, Signatures of the many-body localized regime in two dimensions, Nature Phys 15, 164 (2019).
  * Doggen  _et al._ [2020] E. V. H. Doggen, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Slow many-body delocalization beyond one dimension, Phys. Rev. Lett. 125, 155701 (2020).
  * Kshetrimayum  _et al._ [2020] A. Kshetrimayum, M. Goihl, and J. Eisert, Time evolution of many-body localized systems in two spatial dimensions, Phys. Rev. B 102, 235132 (2020).
  * Nandkishore [2015] R. Nandkishore, Many-body localization proximity effect, Phys. Rev. B 92, 245141 (2015).
  * Luitz  _et al._ [2017] D. J. Luitz, F. Huveneers, and W. De Roeck, How a small quantum bath can thermalize long localized chains, Phys. Rev. Lett. 119, 150602 (2017).
  * Huse  _et al._ [2015] D. A. Huse, R. Nandkishore, F. Pietracaprina, V. Ros, and A. Scardicchio, Localized systems coupled to small baths: From anderson to zeno, Phys. Rev. B 92, 014203 (2015).
  * Wybo  _et al._ [2020] E. Wybo, M. Knap, and F. Pollmann, Entanglement dynamics of a many-body localized system coupled to a bath, Phys. Rev. B 102, 064304 (2020).
  * Hyatt  _et al._ [2017] K. Hyatt, J. R. Garrison, A. C. Potter, and B. Bauer, Many-body localization in the presence of a small bath, Phys. Rev. B 95, 035132 (2017).
  * Goihl  _et al._ [2019] M. Goihl, J. Eisert, and C. Krumnow, Exploration of the stability of many-body localized systems in the presence of a small bath, Phys. Rev. B 99, 195145 (2019).
  * Rubio-Abadal  _et al._ [2019] A. Rubio-Abadal, J.-y. Choi, J. Zeiher, S. Hollerith, J. Rui, I. Bloch, and C. Gross, Many-body delocalization in the presence of a quantum bath, Phys. Rev. X 9, 041014 (2019).
  * Léonard  _et al._ [2020] J. Léonard, M. Rispoli, A. Lukin, R. Schittko, S. Kim, J. Kwan, D. Sels, E. Demler, and M. Greiner, Signatures of bath-induced quantum avalanches in a many-body–localized system, arXiv e-prints , arXiv:2012.15270 (2020), arXiv:2012.15270 [cond-mat.quant-gas] .
  * Krause  _et al._ [2021] U. Krause, T. Pellegrin, P. W. Brouwer, D. A. Abanin, and M. Filippone, Nucleation of ergodicity by a single mobile impurity in supercooled insulators, Phys. Rev. Lett. 126, 030603 (2021).
  * Brighi  _et al._ [2021] P. Brighi, A. A. Michailidis, D. A. Abanin, and M. Serbyn, Propagation of Many-body Localization in an Anderson Insulator, arXiv e-prints , arXiv:2109.07332 (2021), arXiv:2109.07332 [cond-mat.dis-nn] .
  * Verstraete and Cirac [2006] F. Verstraete and J. I. Cirac, Matrix product states represent ground states faithfully, Phys. Rev. B 73, 094423 (2006).
  * Serbyn  _et al._ [2016] M. Serbyn, A. A. Michailidis, D. A. Abanin, and Z. Papić, Power-law entanglement spectrum in many-body localized phases, Phys. Rev. Lett. 117, 160601 (2016).
  * Yu  _et al._ [2017] X. Yu, D. Pekker, and B. K. Clark, Finding matrix product state representations of highly excited eigenstates of many-body localized hamiltonians, Phys. Rev. Lett. 118, 017201 (2017).
  * Khemani  _et al._ [2016] V. Khemani, F. Pollmann, and S. L. Sondhi, Obtaining highly excited eigenstates of many-body localized hamiltonians by the density matrix renormalization group approach, Phys. Rev. Lett. 116, 247204 (2016).
  * Lüschen  _et al._ [2017b] H. P. Lüschen, P. Bordia, S. S. Hodgman, M. Schreiber, S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch, and U. Schneider, Signatures of many-body localization in a controlled open quantum system, Phys. Rev. X 7, 011034 (2017b).
  * Thouless [1972] D. J. Thouless, A relation between the density of states and range of localization for one dimensional random systems, Journal of Physics C: Solid State Physics 5, 77 (1972).
  * Nandy  _et al._ [2021] S. Nandy, F. Evers, and S. Bera, Dephasing in strongly disordered interacting quantum wires, Phys. Rev. B 103, 085105 (2021).
  * Gopalakrishnan  _et al._ [2017] S. Gopalakrishnan, K. R. Islam, and M. Knap, Noise-induced subdiffusion in strongly localized quantum systems, Phys. Rev. Lett. 119, 046601 (2017).
  * Lorenzo  _et al._ [2018] S. Lorenzo, T. Apollaro, G. M. Palma, R. Nandkishore, A. Silva, and J. Marino, Remnants of anderson localization in prethermalization induced by white noise, Phys. Rev. B 98, 054302 (2018).
  * Gefen and Schön [1984] Y. Gefen and G. Schön, Effect of inelastic processes on localization in one dimension, Phys. Rev. B 30, 7323 (1984).
  * Logan and Wolynes [1987] D. E. Logan and P. G. Wolynes, Dephasing and anderson localization in topologically disordered systems, Phys. Rev. B 36, 4135 (1987).
  * Taylor and Scardicchio [2021] S. R. Taylor and A. Scardicchio, Subdiffusion in a one-dimensional anderson insulator with random dephasing: Finite-size scaling, griffiths effects, and possible implications for many-body localization, Phys. Rev. B 103, 184202 (2021).
  * Lezama and Bar Lev [2021] T. L. M. Lezama and Y. Bar Lev, Logarithmic, noise-induced dynamics in the Anderson insulator, arXiv e-prints , arXiv:2108.08861 (2021), arXiv:2108.08861 [cond-mat.dis-nn] .
  * Vidal [2003] G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Phys. Rev. Lett. 91, 147902 (2003).
  * Lukin  _et al._ [2019b] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, and M. Greiner, Probing entanglement in a many-body–localized system, Science 364, 256 (2019b).
  * Kiefer-Emmanouilidis  _et al._ [2020] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Evidence for unbounded growth of the number entropy in many-body localized phases, Phys. Rev. Lett. 124, 243601 (2020).
  * Kiefer-Emmanouilidis  _et al._ [2021] M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Slow delocalization of particles in many-body localized phases, Phys. Rev. B 103, 024203 (2021).
  * Luitz and Lev [2020] D. J. Luitz and Y. B. Lev, Absence of slow particle transport in the many-body localized phase, Phys. Rev. B 102, 100202 (2020).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * Lezama  _et al._ [2017] T. L. M. Lezama, S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, One-particle density matrix occupation spectrum of many-body localized states after a global quench, Phys. Rev. B 96, 060202 (2017).
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, Area laws in a many-body localized state and its implications for topological order, J. Stat. Mech. 2013, P09005 (2013).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), January 2011 Special Issue.