Reviving product states in the disordered Heisenberg chain

Abstract

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.

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 coarse-grained 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²⁰.

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 of 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 that 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 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 high-energy 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: (1) 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; (2) Even given a candidate’s 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 high-fidelity 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.

Results

We focus on the paradigmatic disordered spin-1/2 Heisenberg model on L lattice sites,

$$\hat{H}=-\mathop{\sum }\limits_{j=1}^{L-1}{S}_{j}\cdot {S}_{j+1}+\mathop{\sum }\limits_{j=1}^{L}{h}_{j}{\hat{S}}_{j}^{(z)},$$

(1)

where ${S}_{j}={({\hat{S}}_{j}^{(x)},{\hat{S}}_{j}^{(y)},{\hat{S}}_{j}^{(z)})}^{\top }$ 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³³. 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 that oscillates indefinitely with an amplitude that is lower-bounded by a certified amplitude A_cert. (“Results: Locally oscillating product states”).

Second, we use large-scale tensor network numerics to construct such eigenstates for the disordered Heisenberg chain (“Results: Numerical construction”). 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 (“Results: Main results”).

Lastly, we present theoretical arguments suggesting that large systems may, in fact, host a finite density of locally oscillating excitations (“Results: Multiple localized dynamical oscillations”).

Our results are illustrated in Fig. 1. To keep the discussion concise, we delegate most technical details to the “Methods” section and the Supplementary Information.

Fig. 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 $\langle 2{\hat{S}}_{j}^{x}\rangle$ 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.

Locally oscillating product states

Let us consider two eigenstates $\left|{E}_{1}\right\rangle$ and $\left|{E}_{2}\right\rangle$. Their time-evolved equal superposition

$$\left.\left|{{\Psi }}{(t)}_{\pm }\right.\right\rangle=\frac{1}{\sqrt{2}}\left({{{{{{{{\rm{e}}}}}}}}}^{-{{{{{{{\rm{i}}}}}}}}{E}_{1}t}\left.\left|{E}_{1}\right.\right\rangle \pm {{{{{{{{\rm{e}}}}}}}}}^{-{{{{{{{\rm{i}}}}}}}}{E}_{2}t}\left.\left|{E}_{2}\right.\right\rangle \right)$$

(2)

shows perfect revivals at even multiples of the period τ = π/(E₁ − E₂). Now suppose there is a product state

$$\left| {{\Phi }}{(0)}_{\pm }\right\rangle=\left| {\phi }_{\pm }^{(1)}\right\rangle \otimes \cdots \otimes \left| {\phi }_{\pm }^{(L)}\right\rangle$$

(3)

that approximates $\left|{{{\Psi }}}_{\pm }(0)\right\rangle$ in the sense that its overlap fulfills ${F}_{\pm }^{2}=| \left\langle {{\Psi }}{(0)}_{\pm }| {{\Phi }}{(0)}_{\pm }\right\rangle {| }^{2}\ge 1-\epsilon$ with ϵ small. This implicitly defines the local quantum states $|{\phi }_{\pm }^{(k)}\rangle$. The simple but key observation of our approach is that the time-evolved state $\left|{{\Phi }}{(t)}_{\pm }\right\rangle=\exp (-{{{{{{{\rm{i}}}}}}}}\hat{H}t)\left|{{\Phi }}{(0)}_{\pm }\right\rangle$ will necessarily also show high-fidelity revivals:

$${\left|\left\langle {{\Phi }}{(0)}_{\pm }| {{\Phi }}{(2k\tau )}_{\pm }\right\rangle \right|}^{2}\ge 1-4\epsilon$$

(4)

for any integer k. Moreover, let $j={{{{{{{{\rm{argmin}}}}}}}}}_{k}| \langle {\phi }_{+}^{(k)}| {\phi }_{-}^{(k)}\rangle |$. Then the observable

$$\hat{A}={\mathbb{1}} \otimes \left(\Big| {\phi }_{+}^{(\;j)} \Big\rangle \Big\langle {\phi }_{+}^{(\;j)} \Big| - \Big| {\phi }_{-}^{(\;j)} \Big\rangle \Big\langle {\phi }_{-}^{(\;j)} \Big| \right) \otimes {\mathbb{1}}$$

(5)

is supported on a single site, and its time-dependent expectation value in the state $\left|{{{\Phi }}}_{+}(t)\right\rangle$ oscillates with period τ:

$$\left| \left\langle {{{\Phi }}}_{+}\left| \hat{A}(2k\tau )\right| {{{\Phi }}}_{+}\right\rangle -\left\langle {{{\Phi }}}_{+}\left| \hat{A}((2k+1)\tau )\right| {{{\Phi }}}_{+}\right\rangle \right| \ge {A}_{{{{{{{{\rm{cert.}}}}}}}}}$$

(6)

for any integer k. $\hat{A(t)}$ refers to the Heisenberg picture. The certified amplitude A_cert. is given by

$${A}_{{{{{{{{\rm{cert.}}}}}}}}}=\max \{1-{f}^{2}-2\sqrt{(1-{f}^{2})\epsilon },0\},$$

(7)

where ${f}^{2}=\mathop{\min }\nolimits_{j}| \langle {\phi }_{+}^{(j)}| {\phi }_{-}^{(j)}\rangle {| }^{2}$ measures the minimal local overlap between $|{{\Phi }}{(0)}_{+}\rangle$ and $|{{\Phi }}{(0)}_{-}\rangle$ (assuming that each $| {\phi }_{\pm }^{(j)}\rangle$ is normalized). A detailed proof can be found in “Methods: Certified amplitudes”.

Numerical construction

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 sublattice entanglement entropy (ABABAB. . . -bipartition) is small; see Supplementary Material for more details. Targeting small sublattice entanglement is a heuristic choice motivated by the following considerations: (1) Product states have vanishing sublattice entanglement entropy and therefore any state sufficiently close to a product state should have small sublattice entanglement and (2) even generic translationally invariant matrix-product states (MPS)^69,70, which are commonly considered to be low-entangled, have extensive sublattice entanglement entropies⁷¹. Therefore small sublattice 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 “Methods: Details of our numerical method”).

At sufficiently strong disorder, the eigenstates of $\hat{H}$ feature an area-law entanglement and may be represented faithfully as MPS⁵⁷, whose explicit representation can be determined using the DMRG-X algorithm⁷². The algorithm starts with a “seed” state $\left|{m}_{1}\right\rangle \otimes \cdots \otimes \left|{m}_{L}\right\rangle$, where $| {m}_{j}\rangle \in \{| \uparrow \rangle,| \downarrow \rangle \}$ denote the eigenstates of ${\hat{S}}_{j}^{(z)}$. These seeds are the eigenstates of $\hat{H}$ 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 $\left|E:k\right\rangle$ associated with seeds in domain-wall form

$$\left|{{{{{{{\rm{dw}}}}}}}}:k\right\rangle=\underbrace{\left|\downarrow \right\rangle \otimes \cdots \otimes \left|\downarrow \right\rangle }_{\begin{array}{c}k{{{{{{{\rm{times}}}}}}}}\end{array}}\otimes \left|\uparrow \right\rangle \otimes \cdots \otimes \left|\uparrow \right\rangle .$$

(8)

We then form the superposition of the energy eigenstates resulting from neighboring domain walls,

$$\left| {{{\Psi }}}_{\pm }^{(k)}\right\rangle=\frac{1}{\sqrt{2}}\left(\left| E:k\right\rangle \pm \left|E:k+1\right\rangle \right),$$

(9)

and finally construct their product-state approximation $|{{{\Phi }}}_{\pm }^{(k)}\rangle$. 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 “Methods: Details of our numerical method” for further details). We stress that at this point, it is not clear why the states $|{{{\Psi }}}_{\pm }^{(k)}\rangle$ 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.

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 $\hat{A}$ and that there may exist local operators which oscillate with even higher amplitude.

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 (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 σ²/E² 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.

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 in 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 $\hat{H}$ 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 superpositions 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 $\left|{{{{{{{\rm{dw}}}}}}}}:k\right\rangle$ approach a steady state with a smeared-out interface, a process known as domain-wall melting⁷³. An interface spin pointing away from the z-axis therefore protects against this mechanism.

Since the DMRG-X algorithm outputs the energy eigenstates $\left|E:k\right\rangle$ as MPS, we can compute the expectation values $\langle {{{\Psi }}}_{\pm }^{(k)}| \hat{B(t)}| {{{\Psi }}}_{\pm }^{(k)}\rangle$ of any local operator $\hat{B}$ exactly for arbitrary times t (see “Methods: Details of our numerical method”). This, in turn, allows us to quantitatively estimate the finite-time expectation value $\langle {{{\Phi }}}_{\pm }^{(k)}| \hat{B(t)}| {{{\Phi }}}_{\pm }^{(k)}\rangle$ for any local $\hat{B}$, which is useful since the certified amplitude only provides a lower bound for the oscillations of the specific local observable $\hat{A}$ (yet at infinite times). In Fig. 1, we visualize this for ${\hat{S}}^{(x)}$, which is not strictly identical with the observable $\hat{A}$.

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 ${\hat{S}}_{j}^{(z)}{\hat{S}}_{j+1}^{(z)}$, 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.

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 $\left|{m}_{1}\right\rangle \otimes \cdots \otimes \left|{m}_{L}\right\rangle$ and the energy eigenstates $| {E}_{j}\rangle$ both provide an orthonormal basis of the Hilbert space, there exists a unitary mapping $\hat{U}$ 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 $\hat{U}$ as a local quantum circuit of finite depth and composed of gates that only couple nearest neighbors. At the same time, the Hamiltonian $\hat{H}$, and therefore also the unitary $\hat{U}$, has the states $\left|\downarrow \right\rangle \otimes \cdots \otimes \left|\downarrow \right\rangle$ and $\left|\uparrow \right\rangle \otimes \cdots \otimes \left|\uparrow \right\rangle$ as eigenstates. In the bulk of a large region of spins, all pointing upward or downward, $\hat{U}$ must therefore act like the identity. Quasi-locality immediately implies that $\left|E:k\right\rangle=\hat{U}\left|{{{{{{{\rm{dw}}}}}}}}:k\right\rangle$ only contains a localized, static excitation around the domain-wall interface, see Fig. 3. The superposition $| {{{\Psi }}}_{\pm }^{(k)}\rangle$, which shows perfect revivals, must hence support an operator localized around k whose expectation value oscillates in time, i.e., a dynamical, localized excitation.

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 $\hat{U}$ that maps eigenstates of ${\hat{S}}_{j}^{(z)}$ 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 $\hat{U}$ 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 $\hat{U}$ yields a finite density of localized, dynamical excitations.

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 Fig. 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.

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 ref. ⁷⁵ 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⁸⁰ using the concept of dynamical symmetries⁸¹. 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 we also cannot rule out such a structure. We therefore leave a detailed investigation for future work.

Basic MBL phenomenology has been successfully demonstrated experimentally using ultra-cold Fermions in optical lattices¹⁸ and trapped ions¹⁹. 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 system sizes. Since the preparation of deformed domain walls only requires precise single-site addressing for a 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.

Methods

Certified amplitudes

We derive Eq. (4) and show how to determine the local spin observable $\hat{A}$ that oscillates with the certified amplitude given in Eq. (7). We make use of the general relation

$$| \left\langle {{\Psi }}| {{\Phi }}\right\rangle {| }^{2}=1-D{[\hat{{{\Psi }}},\hat{{{\Phi }}}]}^{2}$$

(10)

between the fidelity and the trace distance D for two pure states, where we use the notation $\hat{{{\Psi }}}=\left.\left|{{\Psi }}\right.\right\rangle \,\,\left.\left\langle {{\Psi }}\right.\right|$. The trace distance fulfills the triangle inequality:

$$ | \left\langle {{{\Phi }}}_{+}(0)| {{{\Phi }}}_{+}({t}_{2k})\right\rangle {| }^{2}\\ \ge 1-{(D[{\hat{{{\Phi }}}}_{+}(0),{\hat{{{\Psi }}}}_{+}(0)]+D[{\hat{{{\Phi }}}}_{+}({t}_{2k}),{\hat{{{\Psi }}}}_{+}(0)])}^{2},$$

(11)

where t_n = nτ. Employing ${\hat{{{\Psi }}}}_{+}({t}_{2k})={\hat{{{\Psi }}}}_{+}(0)$ as well as the fact that the trace distance is invariant under unitary transformations and hence under time-translation, we get $D[{\hat{{{\Phi }}}}_{+}(0),{\hat{{{\Psi }}}}_{+}(0)]=D[{\hat{{{\Phi }}}}_{+}({t}_{2k}),{\hat{{{\Psi }}}}_{+}(0)]=\sqrt{1-{F}_{+}^{2}}$. This yields

$$| \left\langle {{{\Phi }}}_{+}(0)| {{{\Phi }}}_{+}({t}_{2k})\right\rangle {| }^{2}\ge 1-4(1-{F}_{+}^{2})\ge 1-4\epsilon,$$

(12)

where we used the assumption ${F}_{+}^{2}\ge 1-\epsilon$.

We now turn to the operator $\hat{A}$ and its certified amplitude. Let $j={{{{{{{{\rm{argmin}}}}}}}}}_{k}| \langle {\phi }_{+}^{(k)}| {\phi }_{-}^{(k)}\rangle |$ be the site where the local overlap between $\left.\left|{{{\Phi }}}_{-}\right.\right\rangle$ and $\left.\left|{{{\Phi }}}_{+}\right.\right\rangle$ is minimized so that $f=| \langle {\phi }_{+}^{(j)}| {\phi }_{-}^{(j)}\rangle |$. We then define $\hat{A}$ as

$$\hat{A}={\mathbb{1}} \otimes \cdots \otimes \left(\Big| {\phi }_{+}^{(j)} \Big\rangle \Big\langle {\phi }_{+}^{(j)}\Big| -\Big| {\phi }_{-}^{(j)}\Big\rangle \Big\langle {\phi }_{-}^{(j)}\Big| \right) \otimes \cdots \otimes {\mathbb{1}}$$

(13)

The operator-norm of $\hat{A}$ is given by $\parallel \hat{A}\parallel=\sqrt{1-{f}^{2}}$. For any observable $\hat{X}$ and any two density matrices $\hat{\rho }$ and $\hat{\sigma }$ it holds that

$$\begin{array}{r}| {{{{{{{\rm{Tr}}}}}}}}[\hat{X}\hat{\rho }]-{{{{{{{\rm{Tr}}}}}}}}[\hat{X}\hat{\sigma }]| \le \parallel \hat{X}\parallel \,D[\hat{\rho },\hat{\sigma }].\end{array}$$

(14)

Using ${\hat{{{\Psi }}}}_{-}(0)={\hat{{{\Psi }}}}_{+}({t}_{2k+1})$, we therefore find

$$ | {{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{+}({t}_{2k+1})]-{{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Psi }}}}_{-}(0)]| \\ \le \parallel \hat{A}\parallel D[{\hat{{{\Phi }}}}_{+}({t}_{2k+1}),{\hat{{{\Psi }}}}_{+}({t}_{2k+1})]$$

(15)

$$=\sqrt{(1-{f}^{2})(1-{F}_{+}^{2})}$$

(16)

$$\le \sqrt{(1-{f}^{2})\epsilon },$$

(17)

where we used ${F}_{\pm }^{2}\ge 1-\epsilon$. Similarly,

$$| {{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{-}(0)]-{{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Psi }}}}_{-}(0)]| \le \sqrt{(1-{f}^{2})\epsilon }.$$

(18)

The triangle inequality then yields

$$| {{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{+}({t}_{2k+1})]-{{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{-}(0)]| \le 2\sqrt{(1-{f}^{2})\epsilon },$$

(19)

and a similar calculation shows

$$| {{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{+}({t}_{2k})]-{{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{+}(0)]| \le 2\sqrt{(1-{f}^{2})\epsilon }.$$

(20)

Since $f=| \langle {\phi }_{+}^{(j)}| {\phi }_{-}^{(j)}\rangle |$, we further have

$${{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{+}(0)]=1-{f}^{2},\quad {{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{-}(0)]={f}^{2}-1.$$

(21)

In total, we find

$${{{{{{{\rm{Tr}}}}}}}}[\hat{A}{\hat{{{\Phi }}}}_{+}({t}_{2k})] -{{{{{{{\rm{Tr}}}}}}}}\left[\hat{A}{\hat{{{\Phi }}}}_{+}({t}_{2m+1})\right] \\= {{{{{{{\rm{Tr}}}}}}}}\left[\hat{A}{\hat{{{\Phi }}}}_{+}(0)\right]+{{{{{{{\rm{Tr}}}}}}}}\left[\hat{A}({\hat{{{\Phi }}}}_{+}({t}_{2k})-{\hat{{{\Phi }}}}_{+}(0))\right] \\ -{{{{{{{\rm{Tr}}}}}}}}\left[\hat{A}{\hat{{{\Phi }}}}_{-}(0)\right]+{{{{{{{\rm{Tr}}}}}}}}\left[\hat{A}({\hat{{{\Phi }}}}_{-}(0)-{\hat{{{\Phi }}}}_{+}({t}_{2m+1}))\right.\\ \ge (1-{f}^{2})-2\sqrt{(1-{f}^{2})\epsilon }-(\; {f}^{2}-1)-2\sqrt{(1-{f}^{2})\epsilon }\\= 2(1-{f}^{2}-2\sqrt{(1-{f}^{2})\epsilon })$$

(22)

for any $k,m\in {\mathbb{N}}$.

Details of our numerical method

We use a custom implementation of the DMRG-X algorithm, which takes into account the U(1)-symmetry of $\hat{H}$ (the Hamiltonian commutes with the total magnetization in z-direction). 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 MPS representations of excited eigenstates in disordered systems. It starts with an initial MPS called the “seed” (which in our case is a product state on the basis of ${\hat{S}}_{j}^{(z)}$) 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 σ²/E² 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 1, 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 σ² can reach machine precision, while the rescaled σ²/E² can still be larger than 10^−12 if E is smaller than unity. Table 1 contains 6 such states (L = 10, χ = 32).

Fig. 4: Convergence of the DMRG-X algorithm.

Rescaled energy variance σ²/E² along the DMRG-X sweeps for five DMRG-X runs (different initial seeds and disorder realizations) randomly chosen from the full dataset used for Fig. 2. The different panels correspond to system-sizes L = 20, 40, 80, 160 as indicated. The lines show the absolute value ∣σ²/E²âˆ£; missing dots correspond to negative signs (see the main text for details). After each 20 sweeps, the bond dimension χ is increased.

Table 1 Convergence of the DMRG-X algorithm

Due to numerical rounding errors, the energy variance σ² 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 signaling 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 σ² 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 conclusions remain valid for a time, at least of the order of 10⁶ (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.

Fig. 5: Comparison with exact diagonalization.

For each system-size 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 σ²/E² has fallen to at least 10^−12). Given such a state $\left.\left|{{{\Psi }}}_{{{{{{{{\rm{MPS}}}}}}}}}\right.\right\rangle$, we then obtained the closest eigenstate (in terms of overlap) via exact diagonalization $\left.\left|{{{\Psi }}}_{{{{{{{{\rm{ED}}}}}}}}}\right.\right\rangle$ and computed the deviation of the overlap from unity $\delta : \!\!=1-| \left\langle {{{\Psi }}}_{{{{{{{{\rm{MPS}}}}}}}}}| {{{\Psi }}}_{{{{{{{{\rm{ED}}}}}}}}}\right\rangle |$. 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 ({\delta }_{i})$. Then the orange line is given by $\exp (\mu )$ and the size of the upper error bar by $\exp (\mu+s)-\exp (\mu )$. 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.

Finding the product-state approximation

We now explain how to find a product-state approximation to a superposition $\left.\left|{{{\Psi }}}_{\pm }(0)\right.\right\rangle$. Denote by ${\hat{\rho }}_{\pm }^{(j)}$ the reduced density matrix at site j in the state $\left.\left|{{{\Psi }}}_{\pm }(0)\right.\right\rangle$. As with any spin-1/2 density matrix, it may be written as

$${\hat{\rho }}_{\pm }^{(j)}=\frac{1}{2}{\mathbb{1}}+{r}_{\pm }^{(j)}\cdot {S}_{j},$$

(23)

where ${r}_{\pm }^{(j)}$ is the vector that collects the expectation values of the local Pauli operators,

$$\begin{array}{r}{r}_{\pm }^{(j)}=2\left(\begin{array}{r}\left\langle {{{\Psi }}}_{\pm }(0)\right|{\hat{S}}_{j}^{(x)}\left|{{{\Psi }}}_{\pm }(0)\right\rangle \\ \left\langle {{{\Psi }}}_{\pm }(0)\right|{\hat{S}}_{j}^{(y)}\left|{{{\Psi }}}_{\pm }(0)\right\rangle \\ \left\langle {{{\Psi }}}_{\pm }(0)\right|{\hat{S}}_{j}^{(z)}\left|{{{\Psi }}}_{\pm }(0)\right\rangle \end{array}\right).\end{array}$$

(24)

The reduced density matrix is pure if and only if ${r}_{\pm }^{(\;j)}=| | {r}_{\pm }^{(\;j)} | |=1$, and the product state that best approximates each local Pauli expectation value can be obtained by simply normalizing ${r}_{\pm }^{(\; j)}$ to ${\hat{r}}_{\pm }^{(\;j)}={r}_{\pm }^{(\;j)}/{r}_{\pm }^{(\;j)}$. Hence, our product-state approximation is given by $\hat{{{\Phi }}}{(0)}_{\pm }={\otimes }_{j} | {\phi }_{\pm }^{(\; j)}\rangle \langle {\phi }_{\pm }^{(\;j)} |$ with

$$\left| {\phi }_{\pm }^{(j)}\right\rangle \left\langle {\phi }_{\pm }^{(j)}\right|=\frac{1}{2}{\mathbb{1}}+{\hat{r}}_{\pm }^{(j)}\cdot {S}_{j}.$$

(25)

In order to construct to corresponding MPS, we solve the eigenvalue problem of $\frac{1}{2}{\mathbb{1}}+{\hat{r}}_{\pm }^{(j)}\cdot {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]{\sigma }_{j}}$ at site j (with σ_j = ↑, ↓ in our case). The expectation value of an observable $\hat{O}$ supported at lattice site m is then given by

$$\left.\left\langle \psi \right.\right|\hat{O}\left.\left|\psi \right.\right\rangle=\frac{{{{{{{{\rm{Tr}}}}}}}}[(\mathop{\prod }\nolimits_{j=1}^{m-1}{T}_{{\mathbb{1}}}^{[j]}){T}_{O}^{[m]}(\mathop{\prod }\nolimits_{l=m+1}^{L}{T}_{{\mathbb{1}}}^{[l]})]}{{{{{{{{\rm{Tr}}}}}}}}[\mathop{\prod }\nolimits_{j=1}^{L}{T}_{{\mathbb{1}}}^{[\;j]}]},$$

(26)

where the local transfer operator ${T}_{O}^{[j]}$ is defined for any observable $\hat{O}$ supported at site j as

$${T}_{O}^{[j]}=\mathop{\sum}\limits_{{\sigma }_{{j}_{1}},{\sigma }_{{j}_{2}}}{A}^{[j]{\sigma }_{{j}_{1}}}{O}_{{\sigma }_{{j}_{1}}{\sigma }_{{j}_{2}}}{({A}^{[j]{\sigma }_{{j}_{2}}})}^{*}.$$

(27)

We now discuss how to compute a local time-dependent expectation value of a state

$$\left.\left|{{\Psi }}(t)\right.\right\rangle=\mathop{\sum }\limits_{i=1}^{r}{\alpha }_{i}{{{{{{{{\rm{e}}}}}}}}}^{-{{{{{{{\rm{i}}}}}}}}{E}_{i}t}\left.\left|{E}_{i}\right.\right\rangle$$

(28)

in the case where the energy eigenstates $\left.\left|{E}_{i}\right.\right\rangle$ are given as MPS with matrices ${A}_{i}^{[j]{\sigma }_{j}}$ and a bond dimension χ. The state $\left.\left|{{\Psi }}(t)\right.\right\rangle$ can be expressed as an MPS with bond dimension rχ by setting

$${B}^{[j]{\sigma }_{j}}={\oplus }_{i}{A}_{i}^{[j]{\sigma }_{j}}\quad j \, \ne \, m$$

(29)

$${B}^{[m]{\sigma }_{m}}(t)={\oplus }_{i}{\alpha }_{i}{{{{{{{{\rm{e}}}}}}}}}^{-{{{{{{{\rm{i}}}}}}}}{E}_{i}t}{A}_{i}^{[m]{\sigma }_{m}}.$$

(30)

From now on, let ${T}_{O}^{[j]}$ denote the local transfer operators associated with the tensors ${B}^{[j]{\sigma }_{j}}$. Then the time-dependent expectation value takes the form

$$\left.\left\langle {{\Psi }}(t)\right.\right|\hat{O}\left.\left|{{\Psi }}(t)\right.\right\rangle=\frac{{{{{{{{\rm{Tr}}}}}}}}[{T}_{{{{{{{{\rm{left}}}}}}}}}{T}_{O}^{[m]}(t){T}_{{{{{{{{\rm{right}}}}}}}}}]}{{{{{{{{\rm{Tr}}}}}}}}[{T}_{{{{{{{{\rm{left}}}}}}}}}{T}_{{\mathbb{1}}}^{[m]}(t){T}_{{{{{{{{\rm{right}}}}}}}}}]},$$

(31)

where ${T}_{{{{{{{{\rm{left}}}}}}}}}=\mathop{\prod }\nolimits_{j=1}^{m-1}{T}^{[j]}$ and ${T}_{{{{{{{{\rm{right}}}}}}}}}=\mathop{\prod }\nolimits_{l=m+1}^{L}{T}^{[l]}$. Importantly, these left and right transfer operators are independent of t and can be computed once and for all so that all time dependence 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.

Preliminary exact-diagonalization numerics

We performed preliminary small-scale exact-diagonalization numerics targeting small sublattice 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}\rangle$, compute the second Rényi entropy S₂(E_j) of the reduced density matrix associated with every second lattice site (sublattice entanglement).

• 4. Sort the energy eigenstates according to their sublattice entanglement so that S₂(E_j) ≤ S₂(E_k) if j ≤ k.

• 5. For the m eigenstates with the smallest sublattice entanglement and all pairs (E_j, E_k) with j, k = 1, …, m and j < k, construct product-state approximations

  $$\Big| {{{\Phi }}}_{\pm }^{(j,k)}\Big\rangle \, \approx \, \Big| {{{\Psi }}}_{\pm }^{(j,k)}\Big\rangle : \! \!=\frac{1}{\sqrt{2}}\left(\Big| {E}_{j}\Big\rangle \pm \Big| {E}_{k}\Big\rangle \right)$$

  (32)

  and compute the minimum fidelity ${F}^{(i,j)}=\mathop{\min }\nolimits_{\pm }| \langle {{{\Phi }}}_{\pm }^{(j,k)}| {{{\Psi }}}_{\pm }^{(j,k)}\rangle |$, the magnetization profile of $| {{{\Phi }}}_{\pm }^{(j,k)}\rangle$ (local expectation values of the Pauli operators), as well as the associated certified amplitudes. Typically, we chose m = 20.

• 6. Plot fidelities and certified amplitudes and manually inspect the magnetization profile for those states $| {{{\Phi }}}_{\pm }^{(j,k)}\rangle$ with large fidelities and large certified amplitudes. Exemplary data is shown in 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 A_cert. = 0.88 and clearly corresponds to a deformed domain wall.

The outcome of these numerics was a consistent finding of deformed domain walls with large certified amplitudes, which led to the formulation of the DMRG-X-based algorithm directly targeting deformed domain walls.