Reviving product states in the disordered Heisenberg chain

When a generic quantum system is prepared in a simple initial condition, it typically equilibrates toward a state that can be described by a thermal ensemble. A known exception is localized systems that are non-ergodic and do not thermalize; however, local observables are still believed to become stationary. Here we demonstrate that this general picture is incomplete by constructing product states that feature periodic high-fidelity revivals of the full wavefunction and local observables that oscillate indefinitely. The system neither equilibrates nor thermalizes. This is analogous to the phenomenon of weak ergodicity breaking due to many-body scars and challenges aspects of the current phenomenology of many-body localization, such as the logarithmic growth of the entanglement entropy. To support our claim, we combine analytic arguments with large-scale tensor network numerics for the disordered Heisenberg chain. Our results hold for arbitrarily long times in chains of 160 sites up to machine precision.


I. INTRODUCTION
When a large, closed, interacting quantum many-body system is initialized in a simple initial condition, it typically approaches a state that is stationary when only observed with coarse-grained (e.g, local) observables -the system equilibrates [1,2].In addition, the stationary state of the coarsegrained observables is often well-described by statistical (e.g., canonical) ensembles -the system thermalizes [2][3][4][5].While thermalization is a generic phenomenon and aids the theoretical description, it is not inevitable.One of the most intensely debated exceptions is that of many-body localization (MBL), which is realized in interacting quantum models with a sufficiently strong disorder potential [6][7][8][9].Systems exhibiting MBL provide generic examples of non-ergodic systems that fail to thermalize due to a memory of the local initial conditions, yet they are still equilibrating [10,11].Other key features of MBL phases include an unbounded growth of the entanglement during quantum quenches [12][13][14] and peculiar transport properties [15][16][17].There are now a variety of experimental realizations exhibiting signatures of MBL, including, cold atoms [18,19] and photonic systems [20].
The existence of MBL as a stable phase of matter has recently been questioned and it has been suggested that thermalization actually eventually occurs [21][22][23][24][25].However, it is fair to say that a conclusive picture has not yet emerged [26][27][28][29][30][31].A key obstacle is that many studies are based on an exact diagonalization of small systems and might thus not be representative for the behavior in the thermodynamic limit [32,33].Approaching the problem from the perspective of quantum avalanches has been a major recent direction [34][35][36][37][38][39][40][41][42].
Another exception to the rule of equilibration and thermalization was recently discovered: In so-called many-body scarred systems, there exists a relatively small set of initial product states which may show indefinite revivals of the full many-body wavefunction.When the system is initialized in such an initial state, all physical observables (including local * henrik.wilming@itp.uni-hannover.deFIG. 1. Indefinitely oscillating spins.Top: The disordered Heisenberg chain (L = 20) is initialized in a deformed domain wall product state that has an overlap > 0.994 with a superposition of two energy eigenstates.Middle: Under the unitary time evolution, the local spins remain almost uncorrelated and start to oscillate in the region around the domain wall interface.The solid line shows the expectation value of the Pauli-X observable ⟨2 Ŝx j ⟩ at the center spin in the superposed energy eigenstates.The dynamics of the actual product state is within the associated shaded regions due to its large overlap with the superposition of eigenstates.The blue dotted line indicates the certified amplitude, which provides a lower bound for the magnitude of the oscillations in the infinite-time limit.Bottom: Overlay of different snapshots in time of the expectation values of the local spin operators around the domain wall interface, visualized as arrows within their respective Bloch spheres.ones) show periodic oscillations and the system neither thermalizes nor equilibrates [43][44][45][46][47][48][49][50].The revivals of the wavefunction are connected to the existence of a small set of highenergy eigenstates that exhibit atypically low entanglement, dubbed "quantum (many-body) scars".Conversely, if all energy eigenstates are sufficiently entangled, then initial product states generically equilibrate [51][52][53].
In fact, MBL systems also exhibit quantum many-body scarring, and they do so in a most dramatic way: Not just a few, but all high-energy eigenstates have atypically low entanglement, since the entanglement entropy features an area law [54][55][56][57][58][59] instead of a volume law.(This is the generic situation in interacting systems [60][61][62][63][64][65][66][67][68].) To summarize: Many-body scarred systems host a few slightly entangled eigenstates, and these can be sufficient for a complete breakdown of equilibration in certain initial product states.Conversely, all energy eigenstates in MBL systems are low-entangled.This leads to a natural question: Can MBL systems also host initial product states that show high-fidelity revivals of the wavefunction with corresponding local observables that oscillate indefinitely?
If the answer to this question is "yes", then -contrary to current belief -MBL systems do not generally equilibrate from product states and hence also do not thermalize.Moreover, a further hallmark feature of MBL, namely the slow (logarithmic) but unbounded growth of the entanglement entropy, would be violated for these particular initial conditions.
It is, however, unclear how to approach this problem and how to find such initial conditions for a given MBL Hamiltonian.In particular, there are two key difficulties to be overcome: i) The product state might have to be fine-tuned to the details of the Hamiltonian such as the disorder configuration.However, the set of product states is a continuum, so we cannot simply search through all of them.Furthermore, we cannot exploit algebraic structures (such as symmetries) to guide us; ii) Even given a candidate initial state, how could we make sure that it does not equilibrate?In principle, the revival could happen at arbitrary long times, which cannot be accessed analytically or numerically (even for MBL systems).
In this work, we overcome these difficulties and demonstrate that one can find initial product states featuring highfidelity revivals and local observables that oscillate indefinitely.We combine analytical arguments with state-of-the-art tensor network calculations.Importantly, our approach works for arbitrarily long times, and we can treat systems of up to 160 sites with machine precision.

II. RESULTS
We focus on the paradigmatic disordered spin-1/2 Heisenberg model on L lattice sites, where S S S j = ( Ŝ(x) j , Ŝ(y) j , Ŝ(z) j ) ⊤ is the vector of spin-1/2 angular momentum operators at site j.The local magnetic fields h j ∈ [−W, W ] are sampled independently from a uniform distribution; W is the disorder strength.Exact diagonalization of small systems predicts a crossover from an ergodic to an MBL phase around W ∼ 3.5 [33].In the main part of this work, we set W = 8.
First, we show that if we can find two eigenstates whose superposition is well approximated by a product state, then one can construct a local observable which oscillates indefi-nitely with an amplitude that is lower-bounded by a certified amplitude A cert. (Sec.II A).
Secondly, we use large-scale tensor network numerics to construct such eigenstates for the disordered Heisenberg chain (Sec.II B).We present data for systems of up to L = 160 sites and, up to machine precision, provide a rigorous certificate for the indefinite oscillations of a local observable (Sec.II C).
Lastly, we present theoretical arguments suggesting that large systems may in fact host a finite density of locally oscillating excitations (Sec.II D).
Our results are illustrated in Fig. 1.To keep the discussion concise, we delegate most technical details to the Methods (Sec.IV) and the Supplementary Information.

A. Locally oscillating product states
Let us consider two eigenstates |E 1 ⟩ and |E 2 ⟩.Their timeevolved equal superposition shows perfect revivals at even multiples of the period τ = π/(E 1 − E 2 ).Now suppose there is a product state that approximates |Ψ ± (0)⟩ in the sense that its overlap fulfills This implicitly defines the local quantum states |ϕ (k) ± ⟩.The simple but key observation of our approach is that then the time-evolved state |Φ(t) ± ⟩ = exp(−i Ĥt) |Φ(0) ± ⟩ will necessarily also show high-fidelity revivals: for any integer k.Moreover, let j = argmin k |⟨ϕ is supported on a single site and its time-dependent expectation value in the state |Φ + (t)⟩ oscillates with period τ : for any integer k.Â(t) refers to the Heisenberg picture.The certified amplitude A cert. is given by where f 2 = min j |⟨ϕ − ⟩| 2 measures the minimal local overlap between |Φ(0) + ⟩ and |Φ(0) − ⟩ (assuming that each |ϕ (j) ± ⟩ is normalized).A detailed proof can be found in Sec.IV A.
As a next step, we demonstrate how to find pairs of energy eigenstates whose equal superpositions are well approximated by product states.It is reasonable to hypothesize that such states must have a low entanglement with respect to any bipartition.Therefore we performed a structured search on small systems using exact diagonalization and targeting energy eigenstates whose sub-lattice entanglement entropy (ABABAB...-bipartition) is small, see Supplementary Material for more details.Targeting small sub-lattice entanglement is a heuristic choice motivated by the following considerations: i) Product states have vanishing sub-lattice entanglement entropy and therefore any state sufficiently close to a product state should have small sublattice entanglement and ii) even generic translationally invariant matrix-product states (MPS) [69,70], which are commonly considered to be low-entangled, have extensive sub-lattice entanglement entropies [71].Therefore small sub-lattice entanglement heuristically indicates an amount of entanglement that is small even compared to MPS.Our preliminary analysis showed that pairs of energy eigenstates whose equal superpositions are well approximated by product states exist and that one class of them comes in the form of deformed domain walls (see Fig. 1).This knowledge then allows us to devise an efficient tensor-network based algorithm to study large systems, which we now briefly explain (further details may be found in Sec.IV B).

B. Numerical construction
At sufficiently strong disorder, the eigenstates of Ĥ feature an area-law entanglement and may be represented faithfully as MPS [57], whose explicit representation can be determined using the DMRG-X algorithm [72].The algorithm starts with a "seed" state |m 1 ⟩ ⊗ • • • ⊗ |m L ⟩, where |m j ⟩ ∈ {| ↑⟩, | ↓⟩} denote the eigenstates of Ŝ(z) j .These seeds are the eigenstates of Ĥ in the limit of W → ∞.DMRG-X then iteratively determines an (approximate) eigenstate at finite W that is, in a sense, closest to the initial seed.The main numerical control parameter is the so-called bond dimension χ, which we choose so that high-energy eigenstates are obtained up to machine precision.
In our case, we find the energy eigenstates |E : k⟩ associated with seeds in domain-wall form We then form the superposition of the energy eigenstates resulting from neighboring domain-walls, and finally construct their product-state approximation |Φ ± ⟩.This allows us to calculate the certified amplitude A cert. of Eq. (7).All of these operations can be implemented efficiently and accurately in the MPS representation (see Sec. IV B for further details).We stress that at this point it is not clear why the states Ψ (k) ± should be close to product states apart from the fact that we found revolving product states with a similar structure in our small scale exact-diagonalization numerics (see Supplementary Material).Our main results in the next section show that for domain-wall seeds, closeness to a product state is indeed a generic case for sufficiently strong disorder.This in turn immediately implies the non-equilibrating behavior for the associated product states.

C. Main results
In Fig. 2 our aggregated numerical data for the certified amplitude at varying system sizes up to L = 160 and at a disorder strength W = 8 with 100 disorder realizations per system size is depicted (the corresponding fidelities are discussed in Supplementary Note 1).We find median certified amplitudes of the order of 0.7, essentially independent of the system size with decreasing fluctuations as L increases.Moreover, the maximum certified amplitudes for domain-wall states with interface in the middle half of the system (sites k = L/4 to k = 3L/4) slowly increase with system size, with all sampled realizations reaching A cert. > 0.91 for L = 160.The restriction to states with the interface in the middle half of the system excludes states that can be interpreted as being close to single-particle excitations (see below and Supplementary Note 2).We emphasize that the certified amplitude provides a lower bound to the magnitude of the oscillations of Â and that there may exist local operators which oscillate with even higher amplitude.
In a nutshell, Fig. 2 conclusively demonstrates the (generic) existence of initial product states that host high-fidelity revivals and the existence of local, indefinitely-oscillating, observables in a system of up to 160 sites.The overall shape of these product states is of the form of two domain walls separated by a spin pointing roughly in ±x-direction at their interface.Moving away from the interface, the spins still point away from their original ±z-directions, but with decreasing components in the x−y-plane.This is visualized in Fig. 1.As a side remark, we mention that the Hamiltonian Ĥ may also be interpreted as a Hamiltonian of interacting fermions by a Jordan-Wigner transformation.However, in this picture the parity super-selection rule forbids our reviving product states, since they correspond to super-positions of states with different fermion-number parity.
The fact that we find oscillating deformed domain walls is particularly interesting since previous results indicate that the bare domain-wall states |dw : k⟩ approach a steady-state with a smeared-out interface, a process known as domain-wall melting [73].An interface spin pointing away from the z-axis therefore protects against this mechanism.
Since the DMRG-X algorithm outputs the energy eigenstates |E : k⟩ as MPS, we can compute the expectation values ⟨Ψ for any local B, which is useful since the certified amplitude only provides a lower bound for the oscillations of the specific local observable Â (yet at infinite times).In Fig. 1, we visualize this for Ŝ(x) , which is not strictly identical with the observable Â.
Besides the deformed domain-wall states, there exists a second set of reviving product states that exhibit local oscillations.However, these can be interpreted as a single-particle phenomenon arising from Anderson localization and exist irrespective of the strength of the term Ŝ(z) j Ŝ(z) j+1 , see Supplementary Note 2.
In Supplementary Note 3, we further provide numerical data for the certified amplitude and various disorder strength in the range W = 0.5 to W = 8.One can identify a crossover from an ergodic system to a localized system.

D. Multiple localized dynamical oscillations
Our numerical data clearly demonstrates that product states with high-fidelity revivals and locally oscillating observables exist for the disordered Heisenberg model at sufficiently strong disorder.However, our approach only yields states with single dynamical excitations.We now explain our construction qualitatively from a different point of view and argue for the existence of product states with a finite density of such dynamical excitations.
Since the product states |m 1 ⟩ ⊗ • • • ⊗ |m L ⟩ and the energy eigenstates |E j ⟩ both provide an orthonormal basis of the Hilbert space, there exists a unitary mapping Û between the two.The mapping is believed to be quasi-local [9,55,59,74], which implies that it maps local operators to operators whose support is still localized in space with potentially (sub-)exponential tails.As a simplified model for this situation, we may think of Û as a local quantum circuit of finite depth and composed of gates that only couple nearest neighbors.At the same time, the Hamiltonian Ĥ, and therefore also the unitary Û , has the states In the bulk of a large region of spins all pointing upwards or downwards, Û must therefore act like the identity.Quasi-locality immediately implies that |E : k⟩ = Û |dw : k⟩ only contains a localized, static excitation around the domain-wall interface, see Fig. 3.The superposition |Ψ (k) ± ⟩, which shows perfect revivals, must hence support an operator localized around k whose expectation value oscillates in time, i.e., a dynamical, localized excitation.
This discussion suggests that in a large system we may construct multiple domain walls separated by dynamical, localized excitations as long as the size of each domain wall is sufficiently large.A finite density of local dynamical excitations should hence in principle be possible.However, each such excitation doubles the number of energy eigenstates that need to be superposed, and the cost of simulating such situations scales exponentially with the number of excitations.In Supplementary Figure 3 we provide proof-of-principle numerics in a system of size L = 80 with up to three excitations, supporting the general argument described above, see Supplementary Note 4 for more details.to energy eigenstates is (quasi-)local and leaves the states with all spins pointing up or down invariant.Acting on a domain wall, it therefore yields an energy eigenstate with localized static excitation.c) Acting on the two superposed, neighboring domain walls, the unitary Û yields a localized dynamical excitation with perfect revivals.d) If several sufficiently large domain walls are separated by spins pointing in ±x-directions, acting with Û yields a finite density of localized, dynamical excitations.

III. DISCUSSION
Anderson's discovery that a random potential can have strong effects on the transport properties of a free quantum particle was a milestone in condensed matter physics.In the last decade, the fate of Anderson localization in the presence of two-body interactions has received significant attention, and it is believed that generic non-ergodic -so-called many-body localized -systems exist.A key feature of these systems is that simple initial states do not thermalize while local observables still equilibrate.(Some comments on the recent controversy about the existence of the MBL phase can be found in the introduction.) In this work, we provided analytical and numerical arguments that this picture is not correct and that one can construct simple product states that show a complete absence of both thermalization and equilibration.The full many-body wavefunction exhibits high-fidelity revivals and local spin operators oscillate with large amplitudes.We demonstrated this for the prototypical disordered Heisenberg chain via large-scale tensor network numerics for systems of up to L = 160 sites.Our results hold for arbitrary long times up to machine precision.
We also argued that multiple such localized dynamical excitations exist in large systems, giving rise to a picture reminiscent of "Hilbert-space fragmentation" in systems with quantum many-body scars arising from kinematic constraints (see [75] and references therein).Similar results have been found for systems showing so-called "Stark many-body localization", which are translationally invariant systems reproducing much of the MBL phenomenology [76][77][78][79].In this case, local oscillating observables can be proven to exist [80] using the concept of dynamical symmetries [81].In contrast to these disorder-free systems, in our case all of these features depend on the precise disorder realization.Therefore we do not expect a clean, emergent algebraic structure associated with the subspace spanned by states with multiple excitations, but also cannot rule out such a structure.We therefore leave a detailed investigation for future work.
Basic MBL phenomenology has been successfully demon-strated experimentally using ultra-cold Fermions in optical lattices [18] and trapped ions [19].Due to the efficient nature of our algorithm, it is in principle possible to calculate the non-equilibrating product states on the fly given a (quasi-)random disorder realization, even for relatively large systemsizes.Since the preparation of deformed domain-walls only requires precise single-site addressing for few of the spins (with the remaining spins being in large blocks of all up and all down) it should therefore be possible to observe the resulting revivals in present-day or near-future experiments.
Our results were made possible by developing a method to systematically find fine-tuned initial product states.So far, no general and efficient method exists to find product states that resist equilibration and thermalization in general interacting many-body systems.Devising such an approach to study models that are currently believed to be thermalizing is a fruitful future direction.

A. Certified Amplitudes
We derive Eq. ( 4) and show how to determine the local spin observable Â that oscillates with the certified amplitude given in Eq. (7).We make use of the general relation between the fidelity and the trace distance D for two pure states, where we use the notation Ψ = |Ψ ⟩⟨Ψ |.The trace distance fulfills the triangle inequality: where t n = nτ .Employing Ψ+ (t 2k ) = Ψ+ (0) as well as the fact that the trace distance is invariant under unitary transformations and hence under time-translation, we get where we used the assumption F 2 + ≥ 1 − ϵ.We now turn to the operator Â and its certified amplitude.Let j = argmin k |⟨ϕ − ⟩| be the site where the local overlap between |Φ − ⟩ and |Φ + ⟩ is minimized, so that f = |⟨ϕ − ⟩|.We then define Â as The operator-norm of Â is given by ∥ Â∥ = 1 − f 2 .For any observable X and any two density matrices ρ and σ it holds that Using Ψ− (0) = Ψ+ (t 2k+1 ), we therefore find where we used F 2 ± ≥ 1 − ϵ.Similarly, The triangle inequality then yields and a similar calculation shows − ⟩|, we further have In total we find for any k, m ∈ N.

B. Details of our numerical method
We use a custom implementation of the DMRG-X algorithm which takes into account the U (1)-symmetry of Ĥ (the Hamiltonian commutes with the total magnetization in zdirection).The ensuing time evolution of a state which is not an eigenstate of the total magnetization is computed exactly, see below.
The DMRG-X algorithm provides one way to find MPSrepresentations of excited eigenstates in disordered systems.It starts with an initial MPS called the "seed" (which in our case is a product state in the basis of Ŝ(z) j ) and iteratively updates each tensor of the MPS by sweeping through the chain.This is analogous to a ground state calculation, but instead of minimizing the energy in each update step, one picks the eigenstate of the local Hamiltonian that maximizes the overlap with the previous MPS.The bond dimension is increased every 20 sweeps (see Fig. 4); we use values χ = 2, 4, 8, 16, 24, 32 for our main data.The algorithm terminates once the rescaled energy variance σ 2 /E 2 has fallen to at least 10 −12 (E and σ are the bare energy and standard deviation of energy, respectively).As indicated in Fig. 2, we often even find rescaled energy variances below 10 −14 .In Table I, we show how often it is not possible to reach convergence with a maximum bond dimension of χ = 32 in all the calculations resulting in our main result Fig. 2. One should note that for a system of size L = 10, any state can be encoded with a bond dimension χ = 32; however the absolute energy variance σ 2 can reach machine precision, while the rescaled σ 2 /E 2 can still be larger than 10 −12 if E is smaller than unity.Table I contains 6 such states (L = 10, χ = 32).
Due to numerical rounding errors, the energy variance σ 2 may be negative when the calculation has converged to machine precision, even though variances are always positive semi-definite.In such cases, one observes final fluctuations with the same magnitude but differing signs, clearly signalling that the result should be interpreted as zero, see Fig. 4 for examples.
As shown in Supplementary Note 5, a pure state with energy variance σ 2 behaves as an eigenstate for time-scales at least of the order of 1/σ.Hence, the small threshold for the rescaled energy variance of 10 −12 that we use guarantees on its own that all our conclusion remain valid for a time at least of the order of 10 6 (in the chosen units).In Fig. 5 we nevertheless also provide a comparison of the DMRG-X deformed domain-wall states with the closest eigenstates obtained from exact diagonalization for system-sizes up to L = 14, showing excellent agreement in terms of fidelity.

Finding the product state approximation
We now explain how to find a product-state approximation to a super-position |Ψ ± (0) ⟩.Denote by ρ(j) ± the reduced density matrix at site j in the state |Ψ ± (0) ⟩.As any spin-1/2 density matrix, it may be written as where r r r (j) ± is the vector that collects the expectation values of the local Pauli-operators, r r r The reduced density matrix is pure if and only if r  ± .Hence, our product state approximation is given by Φ(0 The individual data-points δi for the various states are shown as light dots.The orange line corresponds to the log-average: let µ and s denote the mean and standard-deviation of log(δi).Then the orange line is given by exp(µ) and the size of the upper error-bar by exp(µ+ s) − exp(µ).We plot the log-average because the sample mean is strongly biased by the comparably few data-points with δi ∼ 10 −10 , whereas the bulk of the data-points lies significantly below 10 −12 for every system-size.
In order to construct to corresponding MPS, we solve the eigenvalue problem of 1 2 1 + r r r (j) ± • S S S j and construct a product state via the local eigenstates associated with the largest eigenvalue.

Long-time simulation using MPS-representations of eigenstates
Let us consider an MPS defined via local tensors A [j]σj at site j (with σ j =↑, ↓ in our case).The expectation value of an observable Ô supported at lattice site m is then given by where the local transfer operator T

[j]
O is defined for any observable Ô supported at site j as We now discuss how to compute a local time-dependent expectation value of a state in the case where the energy eigenstates |E i ⟩ are given as MPS with matrices A [j]σj i and a bond dimension χ.The state |Ψ(t) ⟩ can be expressed as an MPS with bond-dimension rχ by setting From now on let T

[j]
O denote the local transfer operators associated to the tensors B [j]σj .Then the time-dependent expectation value takes the form where . Importantly, these left and right transfer operators are independent of t and can be computed once and for all, so that all timedependence is contained in the local transfer operator T [m] (t).Therefore, it is possible to compute local, time-dependent expectation values at arbitrary times even for large systems.We used this technique to calculate the expectation values in Fig. 1.

C. Preliminary exact-diagonalization numerics
We performed preliminary small-scale exactdiagonalization numerics targeting small sub-lattice entanglement which allowed us to identify domain walls as promising seeds to construct non-equilibrating product states.This procedure consisted of the following steps for systems of sizes L = 8, 10, 12: 1. Sample a disorder realization.
2. Compute all energy eigenstates via exact diagonalization.
3. For each energy eigenstate |E j ⟩, compute the second Rényi entropy S 2 (E j ) of the reduced density matrix associated to every second lattice site (sublattice entanglement).
4. Sort the energy eigenstates accoding to their sublattice entanglement, so that S 2 (E j ) ≤ S 2 (E k ) if j ≤ k.

5.
For the m eigenstates with smallest sublattice entanglement, and all pairs (E j , E k ) with j, k = 1, . . ., m and j < k, construct product state approximations  2 is on the order of machine precision.This entails that all our analytical predictions remain true up to times of the order of the inverse standard deviation of the energy.To see this, note the following standard bound on how much |ψ ⟩ dynamically deviates from being an eigenstate: Since f (0) = 0, we get

FIG. 2 .
FIG.2.Certified amplitudes.Left: Median (light blue) and maximum (light green) of the certified amplitudes that provide a lower bound for the infinite-time oscillations of a local spin observable in a product state corresponding to a deformed domain wall (the median and maximum are taken w.r.t. the different positions of the domain wall; the maximum is restricted to domain walls with an interface in the middle half of the system, i.e., sites L/4 to 3L/4).We present aggregated data for 100 disorder realizations per system-size with disorder strength W = 8 (each point corresponds to one disorder realization).Dark points with error bars show the mean and standard deviation of the associated values.We also plot the median rescaled energy variances σ 2 /E 2 of the eigenstates determined using the DMRG-X algorithm (light red dots) together with their mean and associated variance (dark red).Right: Certified amplitudes (blue) as well as the rescaled energy variances of the two associated eigenstates (red) for all deformed domain walls and a single disorder realization.

FIG. 3 .
FIG. 3. Qualitative picture for dynamical excitations.a) Two neighboring domain walls superpose to a domain wall separated by a spin pointing in ±x-direction.b) At sufficiently strong disorder, the unitary transformation Û that maps eigenstates of Ŝ(z) j

FIG. 4 .
FIG.4.Convergence of the DMRG-X algorithm.Rescaled energy variance σ 2 /E 2 along the DMRG-X sweeps for five DMRG-X runs (different initial seeds and disorder realizations) randomly chosen from the full data-set used for Fig.2.The different panels correspond to sytem-sizes L = 20, 40, 80, 160 as indicated.The lines show the absolute value |σ 2 /E 2 |, missing dots correspond to negative signs (see the main text for details).After each 20 sweeps, the bond dimension χ is increased.

±
|| = 1, and the product state that best approximates each local Pauli expectation value can be obtained by simply normalizing r r r

FIG. 5 .
FIG. 5.Comparison with exact diagonalization.For each systemsize L ∈ {8, 10, 12, 14} we sampled 100 disorder realizations and computed all L non-trivial deformed domain-wall states per disorder realization using DMRG-X in the same way as for our main results (the algorithm terminates once σ 2 /E 2 has fallen to at least 10 −12 ).Given such a state |ΨMPS ⟩, we then obtained the closest eigenstate (in terms of overlap) via exact diagonalization |ΨED ⟩ and computed the deviation of the overlap from unity δ := 1 − | ⟨ΨMPS|ΨED⟩ |.The individual data-points δi for the various states are shown as light dots.The orange line corresponds to the log-average: let µ and s denote the mean and standard-deviation of log(δi).Then the orange line is given by exp(µ) and the size of the upper error-bar by exp(µ+ s) − exp(µ).We plot the log-average because the sample mean is strongly biased by the comparably few data-points with δi ∼ 10 −10 , whereas the bulk of the data-points lies significantly below 10 −12 for every system-size.

FIG. 6 .
FIG.6.Preliminary numerics.Exemplary data for a single disorder realization with W = 8 and a system of L = 10 spins from our preliminary numerics.Left: Fidelities F (i,j) of the first 40 trial states sorted in non-increasing order and their associated certified amplitudes.Right: Magnetization profile in terms of the expectation value of the Pauli-X, Z operators of each lattice site for the third trial state according to the order on the left.The state has fidelity F (i,j) = 0.998, certified amplitude Acert.= 0.88 and clearly corresponds to a deformed domain wall.
SUPPL.FIG.3.Numerical data in support of multiple excitations.Left: For a system of size L = 80 and one disorder realization at W = 8, we determine all eigenstates associated with down-up domain walls at site k1 and up-down domain walls at site k2.For a given distance d = k2 − k1 between their interfaces, we then cut and glue their respective MPS representations at the midpoint to obtain a candidate state with a pair of localized excitations.Blue points show the associated rescaled energy variance (for varying k1 < k2 but fixed d), the orange points are the corresponding median.For sufficiently large d, we obtain energy eigenstates.Right: The same but for triplets of localized excitations separated by d1 = k2 − k1 > 0 and d2 = k3 − k2 > 0. We vary 30 ≤ k2 ≤ 49 and show the corresponding median.Supplementary Note 5: Approximate eigenstates and timescalesOur numerical procedure yields matrix-product states |ψ ⟩ that are approximate eigenstates in the sense that their energy variance Var( |ψ ⟩ , Ĥ) = ⟨ψ | Ĥ2 |ψ ⟩− ⟨ψ | Ĥ |ψ ⟩

TABLE I .
Convergence of the DMRG-X algorithm-For each system size, the table lists the number of initial states (seeds) that have not reached a rescaled energy variance below 10 −12 at a bonddimension χ.