# A route towards engineering many-body localization in real materials

A. Nietner  Helmholtz Zentrum Berlin, 14109 Berlin, Germany  Dahlem Center for
Complex Quantum Systems, Freie Universität Berlin, 14195 Berlin, Germany  A.
Kshetrimayum  Helmholtz Zentrum Berlin, 14109 Berlin, Germany  Dahlem Center
for Complex Quantum Systems, Freie Universität Berlin, 14195 Berlin, Germany
J. Eisert  Dahlem Center for Complex Quantum Systems, Freie Universität
Berlin, 14195 Berlin, Germany  Helmholtz Zentrum Berlin, 14109 Berlin, Germany
B. Lake  Helmholtz Zentrum Berlin, 14109 Berlin, Germany

###### Abstract

The interplay of interactions and disorder in a quantum many body system may
lead to the elusive phenomenon of many body localization (MBL). It has been
observed under precisely controlled conditions in synthetic quantum many-body
systems, but to detect it in actual quantum materials seems challenging. In
this work, we present a path to synthesize real materials that show signatures
of many body localization by mixing different species of materials in the
laboratory. To provide evidence for the functioning of our approach, we
perform a detailed tensor-network based numerical analysis to study the
effects of various doping ratios of the constituting materials. Moreover, in
order to provide guidance to experiments, we investigate different choices of
actual candidate materials. To address the challenge of how to achieve
stability under heating, we study the effect of the electron-phonon coupling,
focusing on effectively one dimensional materials embedded in one, two and
three dimensional lattices. We analyze how this coupling affects the MBL and
provide an intuitive microscopic description of the interplay between the
electronic degrees of freedom and the lattice vibrations. Our work provides a
guideline for the necessary conditions on the properties of the ingredient
materials and, as such, serves as a road map to experimentally synthesizing
real quantum materials exhibiting signatures of MBL.

_Many-body localization_ (MBL) is an intriguing phenomenon characterized by
the absence of thermalization [1] and conservation of information about
initial conditions beyond those of quantum statistical mechanics [2]. Unlike
other forms of ergodicity breaking, MBL is believed to be robust to
Hamiltonian perturbations and therefore, constitute a genuine, stable non-
equilibrium phase of matter. These systems generalize Anderson localization to
interacting systems and are characterized by a number of interesting
properties. They include the absence of particle transport [3], a
characteristic logarithmic growth of entanglement in time [4, 5], highly
excited states obeying an area law for entanglement entropies [6, 7, 8, 9] and
characteristic multi-partite correlations [10, 11], a memory of initial states
in non-equilibrium dynamics [3, 12, 13], and a number of other striking
features. Because of these unique properties violating the eigenstate
thermalization hypothesis [1], MBL as a phase of matter and the mechanisms
leading to it are currently under study to stabilize and explore other
technological and physical applications such as quantum time crystals [14, 15,
16, 17, 18].

Despite its importance and potential applications, the experimental evidence
of MBL has so far basically exclusively been found in synthetic quantum
systems including cold atomic setups [19] and in superconducting circuits
[20]. Both of these setups require precise control and the manipulation of
individual atoms or qubits. At the same time, these approaches are limited to
a relatively small fixed number of particles. This raises the question of
whether properties of MBL systems can be experimentally probed outside such
set-ups in real _quantum materials_. Unlike the cold atomic set-ups, these
systems are actual quantum many body systems in the thermodynamic limit. While
possible signatures of MBL have been found in some real materials, the effect
of electron-phonon coupling remains unclear in these materials [21]. Similar
indirect evidence has been found in driven nuclear spin-chains [22],
disordered super lattices [23] and disordered magnets [24], while conclusive
results are lacking.

Figure 1: Particle imbalance on the radial axis as a measure of localization
for the binary random bond XXZ model with doping strength $\delta=0.3$.
$\theta$ interpolates between the random bond Ising model at
$\theta=0^{\circ}$, the random bond Heisenberg model at $\theta=90^{\circ}$
and the random bond hopping model at $\theta=180^{\circ}$ (cf. the last
paragraph of Section II). The localization of the random Ising model extends
to non-zero values of $\theta$ and transits into the ergodic phase around
$\theta\approx 54^{\circ}$.

In this work, we propose a route towards synthesizing actual materials
exhibiting properties of many-body localization. We suggest that this can be
done by combining two or more related, iso-structural compounds displaying
spin chain magnetism which can be experimentally synthesized in the
laboratory. We investigate in detail how different ratios of the compounds in
the mixing lead to different localization behaviour, discuss the effects of
the choice of the compounds and the impact of the coupling to phonons and
other features and provide a detailed numerical analysis of the different
scenarios. Hence, our work is meant to act as a _blueprint_ of the
experimental synthesis of real materials featuring properties of MBL. We
complement this guide by an intuitive argument and numerical study concerning
the intricate interplay between the electronic structure and the lattice
phonons. In this way, we suggest how signatures of MBL may be robust albeit
the presence of phonons. What is more, this argument may lend a novel
perspective to an open question about the effective decoupling of the spin
degrees of freedom from the phonon background [1, 25, 21].

##  I Theoretical set-up

Figure 2: The top row shows the entanglement entropy $S$ as a function of time
and the bottom row shows the particle imbalance $\mathcal{I}(t)$ at doping
strength $\delta=0.3$ for various disorder strengths $J_{z}=J_{z}^{(0)}$. Sub
figures (a) and (d) correspond to disorder strength $J_{z}=2$, (b), (e) to
disorder strength $J_{z}=6$ and (c), (f) to disorder strength $J_{z}=10$. Blue
shaded regions correspond to the 100 disorder configurations, the disorder
averages are plotted in red.

###  I.1 Review and methods

Many body localization naturally extends Anderson’s idea of order from
disorder from single particle to many body physics: Coupling the local degrees
of freedom to a source of randomness in an adequate way, the ergodic behaviour
of the time evolution of a quantum system is broken. The most studied model
admitting many body localization is given by what we here call the _standard
MBL_ (sMBL) Hamiltonian

$\displaystyle H_{\text{sMBL}}=$
$\displaystyle\sum_{i}\frac{J_{\perp}}{2}\left(s_{i}^{+}s_{i+1}^{-}+s_{i}^{-}s_{i+1}^{+}\right)$
$\displaystyle+$
$\displaystyle\sum_{i}\left(J_{z}s_{i}^{z}s_{i+1}^{z}+Wh_{i}\sigma_{i}^{z}\right)\,,$
(1)

where the real $W$ is the _disorder strength_ and $h_{i}\in[-1,1]$ is a
uniformly and identically distributed random variable.
$s_{i}^{x},s_{i}^{y},s_{i}^{z}$ denote Pauli operators at site $i$, while
$s_{i}^{\pm}:=s_{i}^{x}\pm is_{i}^{y}$. $J_{\perp}$ is the hopping strength
and $J_{z}$, the interaction strength of the magnetic coupling between
neighbouring spins, on a spin chain. At $J_{z}\rightarrow 0$ the model becomes
Anderson localized for any finite value of $W$. While the localizing effects
are destroyed by introducing a finite interaction $J_{z}\neq 0$ at weak
disorder strengths $W$, the system remains localized for sufficiently large
values of $W$. Moreover, it is easy to see that the model becomes diagonal in
the $z$ basis and that the quasi-local constraints of motion become the local
$s_{i}^{z}$ operators in the limit $W\rightarrow\infty$. While
$H_{\text{sMBL}}$ is defined with respect to a continuous disorder field,
discrete disorder fields with only finite levels have also been studied [26,
27, 28].

To proceed, in this work, we discuss effectively one-dimensional models. These
are captured by _matrix product states_ [29, 30] along with the _time evolving
block decimation_ method for our simulations [31, 32], as implemented by the
TeNPy library [33], that allow to keep track of non-equilibrium quantum
systems for not too long times. The use of such _tensor network methods_ [34,
35] comes along with a number of advantages. They are extremely accurate in
their realm of applicability. What is more, entanglement entropies [6] can be
computed with little effort. This feature will be useful in determining
whether a system exhibits properties of MBL or not as the entanglement entropy
in many body localized systems displays an unbounded logarithmic growth in
time [4, 5], starting with a product state that has no entanglement. For a
state vector $\ket{\psi}$, the _entanglement entropy_ with respect to a
subsystem $A$ is given by the von-Neumann entropy $S(\rho):=-\Tr\rho\ln\rho$
of the reduced quantum state on $A$ as

$\displaystyle S_{A}:=S(\rho_{A})$ (2) with the reduced density matrix
$\displaystyle\rho_{A}:=\Tr_{\backslash A}\outerproduct{\psi}{\psi}\,$ (3)

where $\Tr_{\backslash A}$ is the partial trace on the complement of $A$.
Another key property of many body localized systems is the conservation of
classical information that can be quantified by means of the classical
information capacity [2] that in turn can be phrased in terms of a Holevo
quantity. In the case of $H_{\text{sMBL}}$ at strong disorder strength $W$,
this coincides with the information in the computational basis which is
captured by the particle imbalance. The _particle imbalance_ $\mathcal{I}$
measures the amount of the residual, conserved information in the
computational basis and is defined as

$\displaystyle\mathcal{I}(t):=-\frac{1}{n}\sum_{i=1}^{n}(-1)^{i}\expectationvalue{s_{i}^{z}}{\psi(t)}$
(4)

where $\ket{\psi(t)}=\exp(-itH)\ket{\psi_{0}}$ is the time evolved state
vector,
$\ket{\psi_{0}}=\ket{\uparrow,\downarrow,\uparrow,\dots,\uparrow,\downarrow}$
is the Néel state vector and $n$ is the size of the chain. The _quasi-local
constraints of motion_ (qLCOM) [3, 36] of the many body localized fixed point
$H_{\text{sMBL}}(W\rightarrow\infty)$ are the local $s^{z}$ operators.
Therefore, the particle imbalance $\mathcal{I}(t)$ is a constant of motion for
$W\rightarrow\infty$. For finite $W$, it is a measure of the extent to which
local information is preserved. Hence, as the imbalance relaxes to zero for
ergodic systems, it is an indicator for how far from this fixed point the
system is. The combination of the preserved particle imbalance and logarithmic
growth of entanglement highlights core features of MBL as opposed to Anderson
localization: While the classical information which is in our case measured
with respect to the computational basis and is indicated by the particle
imbalance is preserved in many-body localized systems as well as in Anderson
localized systems in time, the quantum correlations as indicated by the
entanglement entropy may still propagate slowly through the system, while they
are preserved in Anderson localized systems. Indeed, while MBL inhibits
particle transport, slow quantum information propagation still takes place
[37, 38].

###  I.2 Doping induced interaction disorder

We now move to discussing the core topic of this work, which concerns ways of
synthesizing and investigating real materials featuring properties of MBL. Our
premises are the following: Firstly, we will consider _discrete disorder_
featuring a finite number of levels only rather than the continuous disorder
strength, $h_{i}$ of $H_{sMBL}$ as specified in Eq. (1). Secondly, instead of
investigating on site disorder, we will consider disorder directly on the
interactions which can be achieved by doping on non-magnetic sites of magnetic
spin chain compounds. While we believe that this model is closer to
experimental realizations, the analysis does carry over qualitatively from on
site disorder since both models have the same MBL fixed point. And thirdly, we
move away from uniformly independently distributed disorder and investigate
disorder distributions that are mimicking the doping process of the real
materials more closely and are hence correlated.

We aim to exploit the inherent randomness in the _doping process_ of two or
more iso-structural spin chain compounds which have related but different
magnetic Hamiltonians as a source of randomness in the system. We consider two
compounds labelled by $0$ and $1$ both described by an effective XXZ
interaction

$\displaystyle H_{\text{XXZ}}^{(j)}=$
$\displaystyle\sum_{i}\frac{J_{\perp}^{(j)}}{2}(s_{i}^{+}s_{i+1}^{-}+s_{i}^{-}s_{i+1}^{+})+\sum_{i}J_{z}^{(j)}s_{i}^{z}s_{i+1}^{z}$
(5)

where $j=0,1$ denotes the respective compound whose hopping strength is
$J_{\perp}^{(j)}$ and interaction strength is $J_{z}^{(j)}$. The real values
of $J_{\perp}^{(j)}$ and $J_{z}^{(j)}$ are fixed by the choice of the
compounds, they cannot be arbitrarily tuned since in reality only a few
compounds will form in the required structure and can be doped into each
other. If we dope these two materials, the resulting Hamiltonian will consist
of strings of couplings of random length stemming from material $0$
alternating with strings of couplings stemming from material $1$. Hence, to a
reasonably good approximation, the resulting Hamiltonian, $H_{\text{doping}}$,
will be of the form

$\displaystyle H_{\text{doping}}=$
$\displaystyle\sum_{i}\frac{J_{\perp}^{c(i)}}{2}(s_{i}^{+}s_{i+1}^{-}+s_{i}^{-}s_{i+1}^{+})+\sum_{i}J_{z}^{c(i)}s_{i}^{z}s_{i+1}^{z}$
(6)

where $J_{\perp}^{c(i)}$ and $J_{z}^{c(i)}$ are distributed randomly according
to the doping strength $\delta$, which gives the ratio of the two materials.
This is reflected in the configuration $c$, where
$c\in\\{0,1\\}^{\times(n-1)}$ is a length $n-1$ bit string with one bit $c(i)$
for each bond $i$ denoting the bond values due to the doping process. A
fraction $\delta$ of these $n-1$ bonds have coupling values $J_{\perp}^{(1)}$
and $J_{z}^{(1)}$ while the remaining fraction $1-\delta$ of $n-1$ bonds have
$J_{\perp}^{(0)}$ and $J_{z}^{(0)}$.

In the laboratory, one has a very precise control over the doping strength. In
fact, this is the parameter with the highest experimental control. Therefore,
it is crucial to take it into the analysis and keep it as a separate parameter
that is able to influence the stability of the MBL state. Accordingly, one has
to sacrifice the independence of the individual bonds in the sampling process.
To see this consider the following situation: if the doping strength is set to
$\delta=0.5$ and the first $n/2$ bits have been assigned the values $c_{i}=1$
for all $i=1$ through $n/2$ then the remaining bits must be assigned the value
$0$. Hence, the individual bits are not distributed independently of each
other. More generally, the disorder can be sampled according to the following
algorithm.

Algorithm 1 Sample doping configuration

0: System size $n\geq 2$ and doping strength $\delta\in[0,1]$

$c\leftarrow(0,\dots,0)\in\\{0,1\\}^{n-1}$

Sample $N=\lceil\delta(n-1)\rceil$ many distinct positions $i\in[n-1]$ and
denote the set by $c^{(1)}$

for $i\in c^{(1)}$ do

$c_{i}\leftarrow 1$

end for

return $c$

Figure 3: In all plots $J_{z}^{(0)}=10$ and $J_{z}^{(1)}=5$ with homogeneous
$J_{\perp}=1$. The top row shows the entanglement entropy $S$ as a function of
time and the bottom row shows the particle imbalance $\mathcal{I}$ as a
function of time for doping strengths $\delta=0.05,\dots,0.9$. The blue shaded
regions are the plots of the individual runs of the $100$ disorder
configurations and the red plot is the ensemble average.

##  II Analysis

Let us first consider a mean field approximation to establish that the
disorder induced by doping indeed has a non-trivial effect. To this end, we
observe that the Hamiltonian $H_{\text{doping}}$ in Eq. (6) has a mean field
approximation with respect to the initial computational basis state vector
$\ket{\psi_{0}}=\ket{j_{0},\dots,j_{n}}$ with
$j_{i}\in\\{\uparrow,\downarrow\\}$

$\displaystyle H_{\text{doping}}^{\text{MF}}=\sum_{i}$
$\displaystyle\frac{J_{\perp}^{c(i)}}{2}(s_{i}^{+}s_{i+1}^{-}+s_{i}^{-}s_{i+1}^{+})+\sum_{i}\sigma(j_{i+1})J_{z}^{c(i)}s_{i}^{z}$

where $\sigma(\uparrow)=1$ and $\sigma(\downarrow)=-1$ quantifies the local
magnetization. In particular, for homogeneous
$J_{\perp}=J^{(0)}_{\perp}=J^{(1)}_{\perp}$ and disordered
$J_{z}$($J^{(0)}_{z}\neq J^{(1)}_{z}$), we observe that the mean field
Hamiltonian resembles the non-interacting Anderson model for localization
lending a first evidence for the non-trivial effect of interaction disorder.

To obtain evidence for the existence of a many body localized phase for the
Hamiltonian $H_{\text{doping}}$ mimicking doped materials, we simulate the
time evolution of the Hamiltonian $H_{\text{doping}}$ in Eq. (6) using matrix
product states and the TEBD algorithm from TeNPy [33] for different system
sizes up to $L=64$ and bond dimensions up to $\chi=1024$ for different
parameters which we discuss below. We compute the particle imbalance and the
entanglement entropy for the setting where we fix the hopping strength
$J_{\perp}^{(i)}=1$ for both $i=0,1$ while the ratio of the interaction
strength was fixed at $J_{z}^{(0)}/J_{z}^{(1)}=2$. The value of the
$J_{z}^{(0)}$ can then thought of as the disorder strength and can be varied
along with the doping strength $\delta$ which parameterizes the fraction of
material labelled $1$ in the compound. While the choice of these parameters
may appear arbitrary, we would like to highlight that the purpose of our work
is to show that the disorder induced in the interaction (bond) of the
Hamiltonian can lead to signatures of MBL. We have also done simulations for
the disorder induced in the hopping strengths of the Hamiltonian but we have
found that this does not play a significant role in the localization (c.f.
Fig. 1). We have also assumed that the interaction strength of one of the
compounds is twice that of the other. This is motivated by the observations
that (1) both interaction strengths should depend on the disorder strength and
(2) both interaction strengths must be sufficiently different to distinguish
from a homogeneous strong Ising anisotropic material. Moreover, as the
microscopic structure of the materials should be similar for the model of the
doping to be justified, we chose a factor of two such that the discrepancy is
not too large. We have found that, for fixed doping strength, the bigger the
difference in the interaction strength of the two materials, the stronger the
localization (data shown in Appendix A.1).

After fixing the doping strength to $\delta=0.3$ and increasing the disorder
strength $J_{z}^{(0)}$, we find a transition from an ergodic phase for
$J_{z}^{(0)}\lesssim 4$ into a localized phase for $J_{z}^{(0)}\gtrsim 4$. As
shown in Fig. 2, for weak disorder strength $J_{z}^{(0)}=2$, the entanglement
entropy grows approximately linearly until it starts oscillating around its
saturation value, while the imbalance immediately drops close to its
equilibrium value $\mathcal{I}_{\infty}\approx 0$. Increasing the disorder
strength $J_{z}^{(0)}=6$ above the transition point we find logarithmic growth
of entanglement typical of MBL while the imbalance saturates around a finite
value $\mathcal{I}_{\infty}\approx 0.8$. Increasing the disorder strength
further $J_{z}^{(0)}=10$ only strengthens the MBL phenomenology: The entropy
grows slower while still asymptotically logarithmic and the saturation value
for the imbalance $\mathcal{I}_{\infty}\approx 0.93$ moves closer to one such
that we conclude to find evidence for MBL.

The analysis of the effect of the doping strength $\delta$ turns out to be
more subtle. In particular, we could not detect a definite doping strength for
which the MBL phenomenology breaks down. However, due to the sparse occurrence
of doped sites in the weak doping regime $\delta<0.1$ and $\delta>0.9$ large
system sizes and a large bond dimension are required to obtain meaningful
results and to avoid finite size or finite entanglement effects as shown in
Fig. 3, we have investigated the dynamical behaviour of the system for
$\delta=0.05,\dots,0.9$ in terms of the entanglement entropy and the particle
imbalance. Starting from $\delta=0.05$ we find that the localizing behaviour
witnessed by the entanglement entropy and the particle imbalance increases
upon increasing $\delta$. However, a peak is reached around $\delta=0.3$ and
the behaviour becomes less pronounced upon increasing $\delta$ beyond
$\delta=0.5$ (data not shown). This can be explained by the dominance of
$J_{z}^{(0)}$: While larger doping strengths $0.6>\delta>0.3$ may induce a
more complex pattern into the doping configuration $c$ leading to stronger
interference effects, the norm of the unperturbed Hamiltonian corresponding to
the terms diagonal in the $z$ basis only decreases as more terms stemming from
material 1 which carry a weight $J_{z}^{(1)}=J_{z}^{0}/2$ are present. This
explains why the localizing effect is far from being symmetric around the
doping strength centre $\delta=0.5$. Because the dominant fraction of terms
carries a weight $J_{z}^{(1)}$ if $\delta>0.5$, the full phenomenology is
comparable to some $\delta<0.5$ at a decreased disorder strength
$\tilde{J}_{z}^{(0)}$. Fig. 3 (a), (b), (f) and (g) (respectively (e) and (j))
display the effects of the entanglement entropy and particle imbalance in the
weak (respectively strong) doping regime.

Interestingly, the particle imbalance does seem to saturate to a finite value
despite the very weak doping strengths $\delta=0.05$ and $\delta=0.1$
(respectively $\delta=0.9$). Similarly, the entanglement entropy does seem to
grow at most logarithmically. However, we note that the particle imbalance has
a smoother transition into the saturated regime compared to the intermediate
doping strengths. Similarly the entanglement entropy seems to be saturated at
a finite value for the weak doping strength $\delta=0.05$ and has a smaller
variance compared to the intermediate doping strengths. This behaviour is very
reminiscent of the undoped XXZ model with a strong Ising anisotropy for finite
system sizes. Within our numerical setting we, therefore, find evidence for
MBL only for doping strengths $0.9\gtrsim\delta\gtrsim.1$. We note that, while
we do not have any evidence for that, the localizing phenomenology at extremal
doping strengths may still reestablish at sufficiently large system sizes.

While doped materials will have different interaction strengths, they will
realistically also have varying hopping strength. Therefore, we have analyzed
the effect of material dependent hopping terms
$J_{\perp}^{(i)}=J_{\perp}\pm\epsilon$ onto the MBL phenomenology. We have
found that the imbalance and the entanglement entropy behave in exactly the
same manner as in the case of homogeneous hopping terms (data not shown). Thus
we conclude that disorder in the hopping has no observable effect on the MBL
phenomenology.

In order to emphasize the importance of a strong anisotropy leaning towards a
very Ising like material, also dealing as a intuition for the material
selection for experimentalists, we have analyzed the Hamiltonian
$H_{\text{doping}}$ (6) with the $J_{z}^{(i)}$’s and $J_{\perp}^{(i)}$’s
varying between random interaction ($J_{\perp}^{(i)}=0$, $J_{z}^{(i)}\neq 0$)
through random Heisenberg ($J_{\perp}^{(i)}=J_{z}^{(i)}$) to random hopping
($J_{\perp}^{(i)}\neq 0$, $J_{z}^{(i)}=0$) Hamiltonians. We set
$J_{z}^{(0)}=2J_{z}^{(1)}=\cos(\theta)$ and
$J_{\perp}^{(0)}=2J_{\perp}^{(1)}=\sin(\theta)$ and fixed $\delta=0.3$. Thus,
$\theta=0^{\circ}$ corresponds to the random Ising point, $\theta=45^{\circ}$
corresponds to the random Heisenberg point and $\theta=90^{\circ}$ corresponds
to the random hopping point. As shown in Fig. 1, we find that strong disorder
alone is not sufficient to lead to localization but that the disorder must be
strong in a local basis such as the computational basis. This is in good
agreement with the theoretical understanding of MBL [39].

##  III Experimental set-up

Figure 4: Numerical results for the compounds of real materials with doping
strength $\delta=0.3$. (a) and (b) shows the results for the compound made of
SrCo2V2O8 and BaCo2V2O8 while (c) and (d) shows the results for the compound
made of CsCoCl3 and CsCoBr3. (a) and (c) show the entanglement entropy. (b)
and (d) show the particle imbalance. Blue shaded regions correspond to the 100
disorder configurations and the red line the corresponding average.

As we have seen in the previous section, in order to obtain disorder induced
MBL it is crucial to combine spin chain compounds that both come with a strong
Ising-like anisotropy while at the same time being sufficiently different in
terms of their $J_{z}$ coupling strength. From an experimental perspective a
further constraint is that the two materials in question must be sufficiently
similar in terms of their crystal structure and space group so that the
structural properties such as lattice parameter and atomic positions vary
smoothly with doping, preserving the spin chain characteristics so that the
resulting material is indeed composed of alternating strings of compound 0 and
1. This already eradicates a number of promising candidate pairs which, while
well matched in terms of their magnetic Hamiltonians, fail to be sufficiently
similar on a microscopic level.

In this section, we are going to compute the same quantities of entanglement
entropy and particle imbalance as in the previous section for two candidate
pairs of compounds. In both cases the two constituent compounds are chosen to
have very similar or identical space groups and their chemical formulas differ
only by one ion. This ion is not the magnetic ion that forms the spin chain
but rather a non-magnetic ion that influences the intra-chain exchange
interactions by slight alterations of the magnetic exchange pathway and the
spacing of the magnetic ions. While both pairs will fail to admit a clear
evidence of MBL we include this analysis as it provides an intuition for the
material selection and the underlying phenomenology.

The first quantum material we investigate is the compound Sr1-δBaδCo2V2O8 made
from the parent compounds SrCo2V2O8 and BaCo2V2O8 which contain screw-chains
of magnetic Co2+ ions that have effective spin-1/2. As shown in Refs. [40, 41,
42], both materials have very similar microscopic structures and correspond to
anisotropic Heisenberg spin-chain anti-ferromagnets. Here, SrCo2V2O8 will
constitute material 0 and can be described by the coupling strengths
$J_{z}^{(0)}=7\text{meV}$ and $J_{\perp}^{(0)}=3.92\text{meV}$ [40]. Likewise
BaCo2V2O8 constitutes material 1 with coupling strengths
$J_{z}^{(1)}=5.8\text{meV}$ and $J_{\perp}^{(1)}=3.074\text{meV}$ [41, 42]. As
shown in Fig. 4 (a) and (b), this combination quickly thermalizes. In
particular, the particle imbalance immediately drops to zero while the
entanglement entropy grows linearly in time. Both clear indications for the
absence of ergodicity breaking. As the ratio between the interaction and the
hopping terms is only of the order $J_{z}^{(0)}/J_{\perp}^{(0)}\approx 1.8$,
this comes at no surprise as the Ising anisotropy is not yet sufficiently
pronounced (c.f. Fig. 2).

We observe the other extreme of non-localized behaviour for the compound
CsCoCl3(1-δ)Br3δ, made of the Ising-Heisenberg anti-ferromagnets CsCoCl3 and
CsCoBr3 [43] where the anisotropic spin chain magnetism also arises from Co2+
ions. The couplings of the effective Hamiltonian are
$J_{z}^{(0)}=13.41\text{meV}$ and $J_{\perp}^{(0)}=1.314\text{meV}$
corresponding to CsCoBr3 and $J_{z}^{(1)}=11.78\text{meV}$ and
$J_{\perp}^{(1)}=1.008\text{meV}$ for CsCoCl3 respectively. Fig. 4 (c) and (d)
show the entanglement entropy and particle imbalance, respectively, of this
material. In contrast to the previous material, we observe a saturation of the
entanglement entropy at a small finite value as well as a saturation of the
particle imbalance. While the saturation of the imbalance seems to indicate
localizing behaviour the saturation of the entanglement entropy rather hints
at a finite size effect (noting that the results are already converged in the
bond dimension). While uncommon for many body localized systems, this finite
size effect is generic for clean systems admitting a strong Ising anisotropy
at finite system sizes. Moreover, in contrast to the results shown in Fig. 2,
we observe only a small variance of the individual configurations around the
their mean in Fig. 4 and the transition into the saturated regime seems less
sharp. We, therefore, suspect the observed saturations of the entanglement
entropy and the particle imbalance shown in the plots do not witness MBL but
resemble finite size effects induced by the strongly pronounced Ising
anisotropy of each of the materials. Moreover, we suspect that the small
variance is due to the small difference of the respective coupling terms
$J_{z}^{(0)}-J_{z}^{(1)}\approx 1.35\text{meV}$ such that the individual
compounds look almost homogeneous and only admit minor differences in their
evolution (c.f. Appendices A.1 and A.2).

We conclude this section with the following observations already elaborated on
in Section II. In order to achieve many body localization in a compound of two
materials it is crucial that both materials are strongly leaning towards an
Ising anisotropy. At the same time, however, both materials need to be
significantly different in the strength of their Ising coupling. We presented
two possible combinations of different materials, each of which only fulfills
one of the two conditions, and observed that neither of them does give rise to
many body localization. Therefore, we believe that the major challenge in the
search for a compound that admits MBL lies in the trade-off between the
similarity of their microscopic structure versus their pronounced difference
in terms of the strength of the (dominating) Ising coupling. While the first
is necessary to ensure that our doping-toy model is a good effective
description, the latter is, as shown above, crucial for the actual MBL
phenomenology.

##  IV Protecting MBL from phonons

It is generally believed that the phonons in a real material constitute a
thermal bath and that the non-zero electron-phonon coupling will lead to
thermalization of the system [1, 44, 45]. For this reason, this aspect
requires particular attention when devising a blueprint for realizing MBL in
real materials. In this section, we will reinvestigate the microscopic
description of the electron-phonon system in order to derive a fresh
perspective on the robustness of the localizing phenomenology. In particular
this analysis lends guidance towards engineering materials in which the
underlying disorder can be used to control and identify the artefacts of the
induced localization.

To begin with let us recap the microscopic model, where without loss of
generality we directly work in the spin rather than in the electron picture.
The Hamiltonian can be recast into the form

$\displaystyle H=H_{\text{spin}}+H_{\text{phonon}}+H_{\text{spin-phonon}}.$
(7)

Here $H_{\text{spin}}$ is the Hamiltonian describing the microscopic spin
degrees of freedom. The microscopic phonon Hamiltonian is described by a
system of harmonic oscillators representing the free bosons

$\displaystyle H_{\text{phonon}}=\sum_{\langle
i,j\rangle}\omega_{i,j}a_{i,j}^{\dagger}a_{i,j}+\sum_{\begin{subarray}{c}\langle\langle
i,j\rangle,\langle
k,l\rangle\rangle\end{subarray}}\kappa_{i,j,k,l}a_{i,j}^{\dagger}a_{k,l},$ (8)

where every bond of the lattice is identified with an oscillator with real
energy $\omega_{i,j}$ and neighbouring oscillators can exchange excitations
via a hopping term $\kappa_{i,j,k,l}$. Note that we use double indices
referring to the bonds in order to align with the notation of the spin
Hamiltonian. We model the interaction between spin and phonon degrees of
freedom as

$\displaystyle H_{\text{spin-phonon}}=\sum_{\langle
i,j\rangle}g_{i,j}\left(a_{i,j}^{\dagger}+a_{i,j}\right)\left(s_{i}^{+}s_{j}^{-}+s_{i}^{-}s_{j}^{+}\right).$
(9)

This is, a hopping of a spin excitation over the bond $\langle i,j\rangle$
couples to the position operator $X_{i,j}\propto(a_{i,j}+a_{i,j}^{\dagger})$
corresponding to the respective oscillator such that the lattice vibrations
modulate the spin transport according to the coupling $g_{i,j}$. For a
homogeneous material $\omega_{i,j}=\omega$ and $\kappa_{i,j,k,l}=\kappa$ such
that $H_{\text{phonon}}$ is diagonal in the momentum basis and the phononic
degrees of freedom can be expected to form a thermal bath. Because the phonon
dynamics is in general faster than the electronic dynamics this gives rise to
an effective heat bath to which the spin or electric degrees of freedom couple
through the interaction term. The corresponding coupling term $g$ has a finite
material dependent value such that the spin degrees of freedom are believed to
thermalize in the long run.

While this analysis gives good results for homogeneous lattices we argue that
it must break down for disordered systems and at sufficiently low temperatures
and sufficiently strong disorder. In addition, the following argument may give
a microscopic explanation to the effective decoupling between electronic and
phononic degrees of freedom as laid out in, e.g., Ref. [21].

###  IV.1 Isolated phonons

We start by observing the effect of the disorder on the microscopic phonon
Hamiltonian. First we note that for a disordered material the oscillator
energies $\omega_{i,j}$ will be not homogeneous, but distributed according to
the bond-disorder. Hence, the isolated phonon Hamiltonian $H_{\text{phonon}}$
of one dimensional systems will be Anderson localized where the localization
length decreases with the underlying disorder strength [46, 47]. The disorder
in the phonon Hamiltonian is also responsible for the inaptness of the
argument made in Ref. [48], which suggest a divergence of the localization
length at low frequencies.

This case is particularly interesting in three (and higher) dimensions. While
there is no generic Andersson localization in three dimensions, there is a
disorder dependent mobility edge. States below some critical energy
$E_{c}(\Delta)$ will be localized while states above this threshold behave in
an ergodic fashion. Moreover, the mobility edge grows with the disorder
strength $\Delta$ where the exact dependence depends on the specific system
[49, 50, 51].

In contrast, two dimensional systems represent a marginal dimension between
the one and three dimensional phenomenology. While all states throughout the
energy spectrum are localized, the localization length of the phonons depends
on their energy in the following way. Below some critical energy scale
$E(\Delta)$ which depends on the disorder strength $\Delta$ the system is
localized exponentially well. Above this critical value, the localization
length depends logarithmically on the energy. Hence, for most experimental
set-ups this will resemble a mobility edge similar to the three dimensional
case fairly well [49, 51, 52].

###  IV.2 One dimensional systems

The Anderson localization of the phononic degrees of freedom comes in nicely
to the study of MBL with phonons. In truly one dimensional systems the
electronic as well as the phononic degrees of freedom are localized. Hence,
assuming a sufficiently short localization length in both systems, translating
to a sufficiently large disorder strength, the interaction between the two
systems is unlikely to induce any transport and we expect the electronic
properties of the full system to be in qualitative agreement with the
properties of the pure electronic system. In particular we believe that the
compound system will remain many body localized.

A rigorous proof of this claim is beyond the scope of this work, but we would
like to point out that the intuitive counter-argument of the phonons forming a
heat bath is no longer valid in the present set-up. Moreover, as both, the
spin and the phonon system are localized by the same disorder, we believe that
the compound system is described by an MBL-like Hamiltonian in the following
way.

The localization of electrons by MBL is believed to be robust to weak local
Hamiltonian perturbations. The Anderson localized phase of the phonons is
thought to become many body localized in the presence of couplings
sufficiently weak compared to the disorder strength and is hence also believed
to be robust. This can be seen as evidence that this stability can extend to
finite coupling strengths that describe the electron-phonon interaction in one
dimensional systems with phonons for disorder strengths sufficiently strong
compared to the electron-phonon coupling. Finally, we would like to point out
that since the phononic degrees of freedom alone do not thermalize, falsifying
our claim would also give new insights into the mechanisms for reversing the
ergodicity breaking introduced by the disorder.

Figure 5: Imbalance $\mathcal{I}$ (a) - (e) and phonon occupation number
$\expectationvalue{n}$ (f) - (j) time series simulated for $8$ spins and $p=1$
phonon mode with cut-of at $k=2^{5}=32$ levels. The spin Hamiltonian
$H_{\text{spin}}$ is the Heisenberg model with random onsite potential with
disorder strength $W=10$. The plots are for 100 disorder realizations, red
curves are the corresponding averages. The phonon energy $\widetilde{\omega}$
and the spin-phonon coupling $g$ are given by $\widetilde{\omega}=10$ and
$g=0.3$ for (a) and (f); $\widetilde{\omega}=10$ and $g=5$ for (b) and (g);
$\widetilde{\omega}=50$ and $g=5$ for (c) and (h); $\widetilde{\omega}=100$
and $g=5$ for (d) and (i) and $\widetilde{\omega}=500$ and $g=5$ for (e) and
(j). While the the system in the weak coupling limit (a) and (f) remains
localized it becomes smeared out as $g$ is increased in (b) and (g). The
system is brought back to a localizing dynamics as the phonon energy
$\widetilde{\omega}$ is increased sufficiently as in (c) and (h) through (e)
and (j). Note that the occupation number in (i) and (j) is re-scaled by
$10^{-2}$ and $10^{-3}$, respectively.

###  IV.3 Three dimensional systems

The analogous chain of thought is less straight-forward for one dimensional
spin systems embedded in a three dimensional material. The localization length
of the phononic degrees of freedom with energies above the mobility edge
$E_{c}(\Delta)$ will diverge. The high energy sector of phonons therefore can
be qualitatively treated similar to a homogeneous lattice, borrowing the
equilibrating behavior. At low energies, sufficiently smaller than
$E_{c}(\Delta)$, however, most phonons are well localized with short
localization length. Hence, with the same argument as in Section IV.2 those
low energy phonons will not induce any transport. Only the high energy phonons
bias the system towards equilibration. To derive a qualitative toy model, we
neglect the low energy well localized phonons with energies
$e<\widetilde{\omega}$ below a critical energy $\widetilde{\omega}$
sufficiently below the mobility edge. Note that the critical energy may not be
defined by the mobility edge but by a qualitatively too large localization
length when compared to the electronic length scales. Hence, we rewrite the
phononic and the interaction Hamiltonians as

$\widetilde{H}_{\text{phonon}}=\sum_{e>\widetilde{\omega}}\widetilde{\omega}_{e}a_{e}^{\dagger}a_{e}$
(10)

and

$\widetilde{H}_{\text{spin-phonon}}=\sum_{e>\widetilde{\omega}}\sum_{\langle
i,j\rangle}\widetilde{g}_{e}\left(a_{e}^{\dagger}+a_{e}\right)\left(s_{i}^{+}s_{j}^{-}+s_{i}^{-}s_{j}^{+}\right).$
(11)

Where we also introduce a toy model for the spin-phonon coupling: Because the
high energy part of the phonon Hamiltonian is approximately diagonal in the
momentum basis a more rigorous treatment would first Fourier transform Eq. (9)
and then truncate the low momentum modes. While this would lead to a more
complicated form of the Hamiltonian if written in real space, the toy model as
defined in Eq. (11) can be expected to have the same equilibrating potential.
In fact, because the latter Hamiltonian couples all modes to all spins in a
geometrically non-local fashion we believe that it even amplifies the
potential equilibration.

Due to the temperature dependent energy distribution we can treat the
corresponding spin-phonon Hamiltonian with an additional high energy cut of
such that only the phonons from within a finite energy band contribute. This
interplay between the temperature dependent distribution of energy and energy
dependent localization length leads to a (diverging) suppression of the
phononic occupation number and hence an effective decoupling of the electronic
from the phononic degrees of freedom at sufficiently low temperatures. For
very low temperatures and large disorder strength we can neglect the phononic
coupling all along. Moreover, as our analysis in Section IV.5 indicates, when
$\widetilde{\omega}$ is sufficiently large compared to $\widetilde{g}$ the
creation of phonons is strongly suppressed. Thus, in this regime the system
will not contain enough high energy phonons to thermalize the spin degrees of
freedom. We thus believe that traces of the MBL phenomenology will endure this
regime. This is phenomenologically in good agreement with experiments made in
Ref. [23].

###  IV.4 Two dimensional systems

Similar as for Anderson localization, the behavior of one dimensional spin
systems embedded in two dimensional lattices is effectively the same as for
those systems embedded in three dimensional lattices. The analysis is thus the
same as in Section IV.3, only that the mobility edge $E_{c}(\Delta)$ is
replaced by the critical energy at which the localization length transits from
a constant to a logarithmic energy dependence. However we do not expect any
major differences in the phenomenology. This is in qualitative agreement with
thin film experiments such as Refs. [53, 21].

###  IV.5 Numerical analysis

The high energy sector of the three dimensional system is unlikely to be
localized. According to our previous analysis, however, the full system may be
non-ergodic due to the localization of the low energy sector. Hence, we expect
insulating effects due to the localization to be visible in experiment in low
to intermediate energy regions similar as in Ref. [23]. We conclude this
section by a numerical small scale simulation via exact diagonalization to
confirm this intuition. To this end we analyze Eq. (12). As we find evidence
for MBL in the doping Hamiltonian as defined in Eq. (6) only for sufficiently
large systems, we use the well established standard MBL Hamiltonian
$H_{\text{spin}}=H_{\text{sMBL}}$ as in Eq. (1). As such, we simulate the
model Hamiltonian

$H_{\text{model}}^{p,k}(W,\widetilde{\omega},\widetilde{g})=H_{\text{sMBL}}(W)+\widetilde{H}_{\text{phonon}}(\widetilde{\omega})+\widetilde{H}_{\text{spin-
phonon}}(\widetilde{g})$ (12)

where we have set the Heisenberg coupling parameters $J_{z}=J_{\perp}=1$. $W$
denotes the disorder strength and the phononic degrees of freedom are
parametrized by $q$ denoting the modes we include and $k$ the truncation
parameter of each bosonic Hilbert space. Due to the numerical resources we
have to trade off between $k$ and $q$. The delocalizing effects of the phonon
bath should be already present in a bath consisting of a single phonon mode,
i.e., $q=1$. Moreover, the dimensionality $k$ of the bath controls the maximal
number of phonons we can encode in that bath. Therefore, it is natural to
investigate the Hamiltonian for the special case $q=1$ and grow $k$ until
convergence in the particle imbalance is reached. While this is a further
simplification it is a good toy model for the system in question which, we
expect, qualitatively inherits the basic properties of the effect in question.
Moreover, while throughout this work we do not consider interactions of the
spin system to environments other than the phonon degrees of freedom, the
effect of the lab temperature, when sufficiently low, may to a crude
approximation be invoked to argue for a cut of in the bosonic degrees of
freedom. We further set the disorder strength $W=10$ in order to find
sufficiently good MBL as we restrict to small systems $n=8$.

For weak couplings $\widetilde{g}\lesssim 0.3$ and arbitrary small oscillator
energy $\widetilde{\omega}$ we find the spin localization to be almost
unaffected by the presence of phonons (c.f. Eq. (5) (a) and (f)). Due to the
small effective phonon energy $\widetilde{\omega}$ the number of phonons may
grow to a non-negligible value. However, the coupling $\widetilde{g}$ is too
weak compared to the robustness of the spin system in order leave a
significant effect on the localization behaviour. Increasing
$\widetilde{g}\approx 5$ we find that the localization of the spin degrees of
freedom vanishes and that the phonon count grows as we increase $k$, even for
small but finite $0<\widetilde{\omega}\lesssim 10$: Evidence for the absence
of MBL in this regime (c.f. Eq. (5) (b) and (g)). Increasing
$\widetilde{\omega}\gtrsim 10$ with $\widetilde{g}=5$ fixed we find that the
spin degrees of freedom are brought back to localization. In particular, the
phonon occupation gets suppressed and the imbalance is stabilized, as
$\widetilde{\omega}\gtrsim 50$ grows sufficiently large (c.f. Eq. (5) (c)-(e)
and (h)-(j)). This confirms our intuition from the previous section, that spin
or electron localization can be made robust against phonons through a
sufficiently high mobility edge.

##  V Conclusion and outlook

In this work, we have studied a novel mechanism to design and synthesize real
materials that are likely to display many body localized phenomena. Moreover,
we have analyzed the effect of the electron-phonon coupling on the
localization phenomenology, to provide evidence that properties of the MBL
phase may be robust under experimentally realistic circumstances. We find
substantial evidence that materials arising out of doping two Heisenberg
magnets with a strong Ising anisotropy indeed gives rise to many body
localization.

In particular, we find that a necessary ingredient is that the two magnets are
at the same time sufficiently Ising anisotropic, and sufficiently different in
terms of their Ising coupling. The difference of Ising couplings of the two
materials together with the pronouncedness of the Ising anisotropy controls
the quality of the MBL phenomenology. While the doping procedure can also lead
to randomness in the hopping term of the resulting material, we find that
randomness in this term does not play a significant role on the localizing
property of the material, provided the hopping terms remain significantly
smaller than the Ising terms of our effective Hamiltonian. When the hopping
terms dominate, our system tends to delocalize. The doping ratio on the other
hand has a similar yet weaker effect on the quality of MBL phenomenology. Our
numerical results indicate that the resulting material may not be many body
localized for combinations with too weak doping ratios.

As a practical example, we have analyzed the likely outcome of doping
different real materials in the lab. The first is the compound Sr1-δBaδCo2V2O8
made from the compounds SrCo2V2O8 and BaCo2V2O8. Our analysis predicts that
the resulting compound is most likely to thermalize. This is because the Ising
anisotropy in the first compound is not sufficiently pronounced. In a second
example, we considered doping between the materials CsCoCl3 and CsCoBr3 [43].
The situation is more complicated here and the results may look localizing for
finite system size. However, since the difference between the Ising coupling
terms of the two materials does not seem to be sufficiently pronounced, the
resulting material CsCoCl3(1-δ)Br3δ is unlikely to display properties of MBL
as evidenced in our numerical analysis.

We have further explored the effect of the electron-phonon coupling on the
many body localized spin phenomenology. To this aim we derived a novel
effective theory describing the spin phonon coupling making small scale
numerical simulations of the global system more accessible. Our main
theoretical insight is that the randomness in the spin Hamiltonian for the
set-up we investigated will be reflected in a random distribution of local
phonon energies and phonon-phonon couplings. Hence, for sufficiently strong
disorders the phonon Hamiltonian will be described by an Anderson localized
Hamiltonian.

Using results for Anderson localization in truly one dimensional systems we
argue that the many body localized phenomenology will persist after coupling
the spin Hamiltonian to the corresponding phonon system for sufficiently
strong disorder.

Three dimensional systmes, where we think of wires of one dimensional spin
chains embedded in a three dimensional lattice, are particularly intricate.
Due to the mobility edge of the three dimensional Anderson localized phonon
system—which is controlled by the disorder strength—there will exist non-local
phonon states. However, due to their large phonon energy their coupling with
spin degrees of freedom will be effectively reduced. We derived an effective
Hamiltonian that represents the spin Hamiltonian coupled to a high energy
phonon and numerically confirmed that a large phonon energy can counter the
effects of a large electron-phonon coupling. We, therefore, believe that for
sufficiently strong disorder, artefacts of the many body localization will
persist after coupling the spin system to the correspondingly disordered
phonon system. In particular, the spin transport should be reduced depending
on the disorder at sufficiently low energies. This is phenomenologically in
good agreement with the results of Ref. [23].

A similar chain of reasoning works in two dimensions, being in good agreement
with observations about thin film experiments made in Ref. [53, 21].

Our work can be understood in two ways. To start with, we provide a road map
for experimentalists for creating materials that carry signs of many body
localized phenomena. In particular we suggest concrete properties for
materials such that the mix of these is believed to display many body
localized phenomena. Our effective theory for the electron-phonon coupled
system further gives a physical intuition that may lend itself towards the
design of experimentally tractable measurement set-ups to confirm or falsify
the claims made.

In a different reading, we re-analyze the effect of the electron-phonon
coupling present in realistic materials on MBL. Contrary to the common belief,
however, we arrive at the conclusion that meaningful artefacts of many body
localization are actually possible in real materials. In this sense, our work
can be understood as pointing out concrete properties and assumptions that
remain to be analyzed in more detail in order to arrive at a final answer as
to whether many body localization or artifacts thereof can survive when
transiting from the theoretical description to an actual realization of such
systems in real materials.

Either way, our work opens the possibility and creates a clear map for future
research to settle the question of whether real materials can exist carrying
signs of many body localization in a more rigorous manner. It would also be
interesting to consider situations where we mix more than two ingredient
compounds giving rise to higher levels of discrete disorder. Another
possibility is to consider spin chains of higher spins in order to increase
the range of candidate materials that can give rise to real materials with
properties of MBL. It is the hope that the present work can bring MBL in real
quantum materials substantially closer to reality.

## Acknowledgements

We would like to acknowledge discussions with Steven Thomson, Christian
Bertoni, Marcel Goihl, Curt von Keyserlingk, Christian Krumnow and Bohan Lu.
A. N. is very thankful to the local IT team and in particular to Jörg Behrmann
and Jens Dreger for generous support regarding the numerical implementation
and maintenance of the computational resources. Our tensor network simulations
were performed using the TeNPy library [33]. We thank the DFG (CRC 183 and FOR
2724) for funding. This work also touches upon the themes of BMBF (MUNIQC-
Atoms) in the tensor network method development.

Figure 6: Entanglement entropy (top row) and particle imbalance (bottom row)
for different interaction gaps $J_{z}^{0}-J_{z}^{1}$ at doping strength
$\delta=0.3$ and $J_{\perp}^{i}=1$. In all plots $J_{z}^{0}=10$ is fixed and
$J_{z}^{1}$ is varied as $J_{z}^{1}=1$ in (a), (e), $J_{z}^{1}=3$ in (b), (f),
$J_{z}^{1}=7$ in (c), (g), and $J_{z}^{1}=9$ in (d), (h).

##  Appendix A Supplemental material to the numerical investigation

###  A.1 On the role of the interaction strength gap

In this appendix, we investigate the role of the interaction strength gap
$\Delta=J_{z}^{0}-J_{z}^{1}$. Figure 6 shows the data for fixed doping
strength $\delta=0.3$, $J_{\perp}^{i}=1$ and $J_{z}^{0}=10$ for varying
$J_{z}^{1}=1,3,7,9$. The variance of the individual plots becomes vanishing as
$\Delta$ becomes small as in Fig. 6 (d) and (h) ($J_{z}^{1}=0.9$). This
indicates that the corresponding system behaves approximately like the
homogeneous, undoped, material as analyzed in Appendix A.2.

Figures 6 (b), (c), (f) and (g) resemble the corresponding plots for
$J_{z}^{1}=5$ as in the main text and can be seen as evidence for MBL. As
$J_{z}^{1}$ becomes as small as $1$ in Figs. 6 (a) and (e) the effect becomes
less pronounced. We expect that this is due to the proximity of the hopping
strength and the interaction strength of the second material.

###  A.2 Finite size scaling for strong Ising anisotropy

In this appendix, we analyze the homogenous, i.e., undoped XXZ Hamiltonian at
strong Ising anisotropy. As shown in Fig. 7 we find that the the entanglement
entropy as well as the particle imbalance saturates at any fixed system size
to a finite value. This value, however, increases, respectively decreases, as
the system size is increased. In particular, we do not expect to find a
saturating particle imbalance in the thermodynamic limit. This agrees with the
fact that the homogenous XXZ model at any finite anisotropy is ergodic.

Figure 7: Finite size scaling for the entanglement entropy (a) and the
imbalance (b) in the homogeneous strong Ising anisotropic XXZ Hamiltonian.

## References

  * Abanin et al. [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Nico-Katz et al. [2022] A. Nico-Katz, A. Bayat, and S. Bose, Phys. Rev. B 105, 205133 (2022).
  * Huse et al. [2014] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
  * Bardarson et al. [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * Žnidarič et al. [2008] M. Žnidarič, T. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).
  * Eisert et al. [2010] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, J. Stat. Mech. P09005 (2013).
  * Luitz et al. [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).
  * Friesdorf et al. [2015a] M. Friesdorf, A. H. Werner, W. Brown, V. B. Scholz, and J. Eisert, Phys. Rev. Lett. 114, 170505 (2015a).
  * Goold et al. [2015] J. Goold, C. Gogolin, S. R. Clark, J. Eisert, A. Scardicchio, and A. Silva, Phys. Rev. B 92, 180202 (2015).
  * De Tomasi et al. [2017] G. De Tomasi, S. Bera, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 118, 016804 (2017).
  * Eisert et al. [2015] J. Eisert, M. Friesdorf, and C. Gogolin, Nature Phys. 11, 124 (2015).
  * Polkovnikov et al. [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Wilczek [2012] F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * Else et al. [2016] D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Yao et al. [2017] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Zhang et al. [2017] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath, et al., Nature 543, 217 (2017).
  * Choi et al. [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, et al., Nature 543, 221 (2017).
  * 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, Science 349, 842 (2015).
  * Roushan et al. [2017] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, et al., Science 358, 1175 (2017).
  * Ovadia et al. [2015] M. Ovadia, D. Kalok, I. Tamir, S. Mitra, B. Sacépé, and D. Shahar, Scientific Rep. 5, 13503 (2015).
  * Wei et al. [2018] K. X. Wei, C. Ramanathan, and P. Cappellaro, Phys. Rev. Lett. 120, 070501 (2018).
  * Nguyen et al. [2021] T. Nguyen, N. Andrejevic, H. C. Po, Y. Tsurimaki, N. C. Drucker, A. Alatas, E. E. Alp, B. M. Leu, A. Cunsolo, Y. Q. Cai, et al., Nano Lett. 21, 7419 (2021).
  * Silevitch et al. [2019] D. M. Silevitch, C. Tang, G. Aeppli, and T. F. Rosenbaum, Nature Commun. 10, 4001 (2019).
  * Altshuler et al. [2009] B. L. Altshuler, V. E. Kravtsov, I. V. Lerner, and I. L. Aleiner, Phys. Rev. Lett. 102, 176803 (2009).
  * Andraschko et al. [2014] F. Andraschko, T. Enss, and J. Sirker, Phys. Rev. Lett. 113, 217201 (2014).
  * Kshetrimayum et al. [2020] A. Kshetrimayum, M. Goihl, and J. Eisert, Phys. Rev. B 102, 235132 (2020).
  * Kshetrimayum et al. [2021] A. Kshetrimayum, M. Goihl, D. M. Kennes, and J. Eisert, Phys. Rev. B 103, 224205 (2021).
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Perez-Garcia et al. [2007] D. Perez-Garcia, F. Verstraete, M. Wolf, and J. Cirac, Quant. Inf. Comp. 7, 401 (2007).
  * Vidal [2004] G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
  * Daley et al. [2004] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. p. P04005 (2004).
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, SciPost Phys. p. 5 (2018).
  * Orús [2014] R. Orús, Ann. Phys. 349, 117 (2014).
  * Verstraete et al. [2008] F. Verstraete, J. I. Cirac, and V. Murg, Adv. Phys. 57, 143 (2008).
  * Serbyn et al. [2013] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013).
  * Friesdorf et al. [2015b] M. Friesdorf, A. H. Werner, M., J. Eisert, and W. Brown, New J. Phys. 17, 113054 (2015b).
  * Bañuls et al. [2017] M. C. Bañuls, N. Y. Yao, S. Choi, M. D. Lukin, and J. I. Cirac, Phys. Rev. B 96, 174201 (2017).
  * Imbrie [2016] J. Z. Imbrie, Phys. Rev. Lett. 117, 027201 (2016).
  * Bera et al. [2017] A. K. Bera, B. Lake, F. H. L. Essler, L. Vanderstraeten, C. Hubig, U. Schollwöck, A. Schneidewind, and D. L. Quintero-Castro, Phys. Rev. B 96, 16 (2017).
  * Grenier et al. [2015a] B. Grenier, S. Petit, V. Simonet, E. Canévet, L.-P. Regnault, S. Raymond, B. Canals, C. Berthier, and P. Lejay, Phys. Rev. Lett. 114, 017201 (2015a).
  * Grenier et al. [2015b] B. Grenier, S. Petit, V. Simonet, E. Canévet, L.-P. Regnault, S. Raymond, B. Canals, C. Berthier, and P. Lejay, Phys. Rev. Lett. 115, 119902 (2015b).
  * Nagler et al. [1983] S. E. Nagler, W. J. L. Buyers, R. L. Armstrong, and B. Briat, Phys. Rev. B 27, 1784 (1983).
  * Agarwal et al. [2017] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan, D. A. Huse, and M. Knap, Ann. Phys. 529, 1600326 (2017).
  * Wu et al. [2019] L.-N. Wu, A. Schnell, G. D. Tomasi, M. Heyl, and A. Eckardt, New J. Phys. 21, 063026 (2019).
  * Anderson [1958] P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Abrahams [2010] E. Abrahams, _50 Years of Anderson Localization_ (World Scientific, 2010).
  * Gurarie and Chalker [2003] V. Gurarie and J. T. Chalker, Phys. Rev. B 68, 134207 (2003).
  * Abrahams et al. [1979] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
  * Kondov et al. [2011] S. S. Kondov, W. R. McGehee, J. J. Zirbel, and B. DeMarco, Science 334, 66 (2011).
  * Licciardello and Thouless [1978] D. C. Licciardello and D. J. Thouless, J. Phys. C 11, 925 (1978).
  * White et al. [2020] D. H. White, T. A. Haase, D. J. Brown, M. D. Hoogerland, M. S. Najafabadi, J. L. Helm, C. Gies, D. Schumayer, and D. A. W. Hutchinson, Nature Comm. 11, 4942 (2020).
  * Ovadia et al. [2009] M. Ovadia, B. Sacépé, and D. Shahar, Phys. Rev. Lett. 102, 176802 (2009).