# Many-body (de)localization in large quantum chains

## Abstract

We theoretically study the quench dynamics for an isolated Heisenberg spin chain with a random on-site magnetic field, which is one of the paradigmatic models of a many-body localization transition. We use the time-dependent variational principle as applied to matrix product states, which allows us to controllably study chains of a length up to spins, i.e., much larger than that can be treated via exact diagonalization. For the analysis of the data, three complementary approaches are used: (i) determination of the exponent which characterizes the power-law decay of the antiferromagnetic imbalance with time; (ii) similar determination of the exponent which characterizes the decay of a Schmidt gap in the entanglement spectrum, (iii) machine learning with the use, as an input, of the time dependence of the spin densities in the whole chain. We find that the consideration of the larger system sizes substantially increases the estimate for the critical disorder that separates the ergodic and many-body localized regimes, compared to the values of in the literature. On the ergodic side of the transition, there is a broad interval of the strength of disorder with slow subdiffusive transport. In this regime, the exponents and increase, with increasing , for relatively small but saturate for , indicating that these slow power laws survive in the thermodynamic limit. From a technical perspective, we develop an adaptation of the “learning by confusion” machine learning approach that can determine .

## I Introduction

Many-body localization (MBL) refers to the localization of particles by disorder in the presence of interactions at non-zero energy density (for recent reviews, see Nandkishore and Huse (2015); Altman and Vosk (2015); Abanin and Papić (2017)), as opposed to the conventional Anderson localization Evers and Mirlin (2008) which describes non-interacting particles in the presence of disorder. MBL is of fundamental interest for understanding metal-insulator transitions and disordered superconductors. From a numerical perspective, MBL is a notoriously difficult problem to describe because of its many-body nature, sensitivity to finite-size effects and the requirement for ensemble averaging over many realizations of disorder. On a conceptual level, the study of MBL ties to the thermalization of quantum systems and the bridge between microscopic dynamics and quantum statistics Rigol et al. (2008).

The study of MBL did not begin in earnest until the development of landmark theories Gornyi et al. (2005); Basko et al. (2006) predicting a temperature-driven transition to a localized phase, now known as an MBL phase. At the same time, advances in the increase of computational power and the development of new algorithms based on tensor networks like matrix product states (MPS) have dramatically accelerated the numerical study of disordered interacting systems, albeit thus far restricted to ground states, short times, boundary-driven systems or very strong disorder Žnidarič et al. (2016); Khemani et al. (2016); Yu et al. (2017); Wahl et al. (2017); Doggen et al. (2017). Experimentally, MBL has been reported in ultracold atoms Schreiber et al. (2015); Choi et al. (2016), trapped ions Smith et al. (2016), and dipolar spins in diamond Kucsko et al. (). Analytically, the existence of an MBL region in the phase diagram of a disordered quantum spin chain has been proven based on an assumption about energy level repulsion Imbrie (2016).

In this work, we investigate a Heisenberg chain with isotropic interactions and a random magnetic field, which has become one of the paradigmatic models for the investigation of the MBL-related physics. In particular, Ref. Pal and Huse (2010) considered system sizes from to and found crossing points in level statistics and relaxation of spin modulations plotted as a function of disorder (in notations of the present work) for various system sizes. By analogy with the scaling analysis of the Anderson transition in non-interacting systems, these crossing points serve as an indication of the MBL transition. However, a rather strong drift of the crossing points was found, from to , as increased from 8 to 16, making it difficult to locate the transition. In a later work Luitz et al. (2015), systems of larger sizes, up to , were investigated via exact diagonalization (for the current state-of-the-art, see Ref. Pietracaprina et al. (2018)), and the estimate for the transition point was obtained, consistent with the upper bound found in Ref. Berkelbach and Reichman (2010). On the other hand, Ref. Devakul and Singh (2015) used a numerical linked-cluster expansion for the entanglement entropy in the thermodynamic limit and obtained a lower bound which appears to be in conflict with the exact-diagonalization results quoted above.

All in all, the current status of the numerically obtained results is still controversial; simulations of small systems can easily miss physics relevant at longer length scales even if ameliorated by finite-size scaling approaches. Most importantly, the strong dependence of the apparent transition point on the system size calls for a detailed numerical study of physical observables characterizing the MBL transition in systems of size that is well beyond the reach of exact diagonalization. To achieve this goal, we consider the MPS-based approach, which allows us to controllably study the long-time quench dynamics around the MBL transition in spin chains of length up to . We also apply machine-learning techniques to further analyze the data obtained from MPS.

Clearly, not only the position of the MBL transition is of interest but also a more detailed insight into the physics of phases around it. Previous works on one-dimensional (1D) systems indicated that a part of the delocalized phase adjacent to the transition is characterized by a slow, subdiffusive dynamics Serbyn et al. (2013); Bar Lev et al. (2015); Agarwal et al. (2015, 2016); Luitz et al. (2016). However, it was also found Bera et al. (2017), within exact-diagonalization studies of systems with size up to , that an apparent exponent characterizing this subdiffusive phase changes substantially with the system size. This poses a question of whether the subdiffusive behavior is a genuine property of the system in the limit, or merely a transient feature which characterizes the relatively small systems. To answer this question, it is important to numerically study substantially larger systems. This is another motivation for the present investigation of the MBL physics in a Heisenberg chain within the MPS-based approach.

## Ii Model and method

### ii.1 Disordered Heisenberg XXZ chain

We consider the Heisenberg XXZ chain with an on-site random field on a lattice of length with open boundary conditions, as described by the Hamiltonian

(1) |

Here , , and are the standard spin- Pauli operators corresponding to the site , and the on-site field takes random values according to a uniform distribution . In the following, we set as a choice of units, and put , i.e., we consider the isotropic Heisenberg chain, unless stated otherwise. Using a Jordan-Wigner transformation, this model can be mapped to the model of nearest-neighbor interacting hard-core bosons with an on-site disordered potential. We consider the zero spin sector , which corresponds to half-filling in the particle picture. The Anderson model of non-interacting particles in a disordered potential is retrieved for . The model (1) was used to discuss the MBL transition and the subdiffusive behavior in the delocalized phase, see the references above.

### ii.2 Time-dependent variational principle

To study the model (1), we employ a recently developed numerical method for describing the time evolution of 1D lattice systems, which is based on the Dirac-Frenkel time-dependent variational principle (TDVP) Dirac (1930) as applied to MPS Haegeman et al. (2011, 2016). Contrary to traditional time-evolution algorithms, such as the time-dependent density matrix renormalization group (t-DMRG) Daley et al. (2004); White and Feiguin (2004) or time-evolving block decimation (TEBD) Vidal (2003), the TDVP does not rely on a Suzuki-Trotter decomposition of local Hamiltonian terms. Instead, the time-dependent wave function is given by

(2) |

where projects the time-evolved wave function back onto the variational MPS manifold, with a dimension typically much smaller than the dimension of the complete Hilbert space .

From the perspective of accuracy of MPS simulations, the worst-case scenario for a local Hamiltonian is realized when the von Neumann entanglement entropy for a bipartition of the chain into two parts grows linearly Calabrese and Cardy (2005), which is known as volume-law entanglement and should be contrasted to the area-law entanglement characteristic of localized systems Laflorencie (2016). It is well known that in such a worst-case scenario, the required bond dimension (which controls the dimension of the variational manifold) grows exponentially in time if it is required that the truncation error is kept below some finite value. This provides a fundamental limit to the maximum time reached using MPS simulations under the condition that Schollwöck (2011). Within the traditional DMRG framework, recent developments have extended the maximum possible times significantly Karrasch et al. (2013), although the fundamental limit remains. On the other hand, under certain conditions it is expected Leviatan et al. (2017); Kloss et al. (2018); Berta et al. (2017) that the TDVP can provide reasonable estimates for the transport properties despite a potentially large truncation error. A key advantage of the method is that unlike t-DMRG and TEBD, it conserves the total energy by construction.

In the model (1), an MBL regime characterized by area-law entanglement has been identified for strong disorder (large ). On the other hand, for weak disorder (small ) an ergodic regime has been reported, described by the eigenstate thermalization hypothesis (ETH) Deutsch (1991); Srednicki (1994) and exhibiting volume-law entanglement. In the present work, we focus on the range of sufficiently strong disorder , which includes the whole MBL phase as well as a part of the delocalized phase. As we have verified (see detailed explanations below), the MPS approach works reliably up to long times considered in our study not only in the MBL phase but also in the delocalized phase with since the corresponding dynamics is very slow.

## Iii Numerical results (TDVP)

In this section, we will detail results obtained using the direct analysis of the TDVP data for the quench dynamics in disordered Heisenberg chains (1) of lengths and 100 with disorder strength ranging from to .

### iii.1 Imbalance

We follow the time dynamics of an initial unentangled (product) Néel state , where we consider many different realizations of disorder. Whether the system is in the delocalized or localized phase can be quantified using the imbalance:

(3) |

a quantity that measures how much of the initial antiferromagnetic order remains at time Luitz and Bar Lev (2017). By definition, , and for a fully ergodic, thermalized state, at long times on average, which means that the system completely loses memory of the initial antiferromagnetic order. The disorder average is denoted as , where labels the different realizations.

The imbalance is appealing from several perspectives. It is a global characteristic of the system which signifies its degree of localization, but at the same time it is computed from purely local quantities (average spins). In view of this locality, the imbalance can be readily measured in an experimental setting using cold atoms Schreiber et al. (2015). We can also compare directly the localization properties of an interacting system to those of non-interacting Anderson insulators. The non-interacting value is obtained by computing the time evolution of the Néel state through the exact diagonalization of (1) with (cf. Ref. Berkelbach and Reichman (2010)). In Fig. 1, we show the dynamics of imbalance for both the interacting and non-interacting cases using the parameters , .

For a thermalized system, a power-law decay is expected Luitz and Bar Lev (2017). In fact, within the Boltzmann equation, one would get an exponential decay with time, since the imbalance corresponds to a mode with a large wave vector (namely, ) describing, in contrast to the total spin, a non-conserved quantity. However, taking into account the coupling of this mode to the low- diffusive (or subdiffusive) mode associated with the spin density spreading suggests a power-law decay (by analogy with the long-time tails related to the return probability that are found for other observables). Indeed, previous numerical results indicated that the exponent of the imbalance decay is qualitatively similar to the subdiffusion exponent found from the return probability, the mean-square displacement of a spin excitation, and the low-frequency dependence of the conductivity (see Ref. Luitz and Bar Lev (2017) and references therein).

In Fig. 2, we show the exponent computed numerically using a least-squares fitting algorithm for various values of and from the imbalance at sufficiently long times . We choose as the weakest disorder to be considered. This is based on the requirement that the value of the bond dimension used in our computations is sufficient to ensure that the value of is insensitive to a further increase of (see Appendix A). As a further check, we compute using a completely independent implementation of the time-evolving block decimation (TEBD) algorithm for , , , finding agreement within error bars. In addition, for a few values of and for the system size we follow the time evolution until using and find a good agreement with obtained in the window (not shown), as well as checking that for bond dimension up to , is insensitive to increasing the time window up to for and (see Appendix A). The number of different disorder realizations for each choice of and is typically of the order of .

For small systems that can be treated by exact diagonalization (such as ), we find a good agreement with previous results Luitz et al. (2016). Specifically, the value of , which continuously decreases with increasing , vanishes within error bars at . However, an increase of the system size from to leads to a substantial increase of , consistent with a trend observed for relatively small systems in Ref. Bera et al. (2017). As a result, the value of disorder for which vanishes also increases
^{1}

We thus conclude that, after the initial increase for , the exponent saturates, i.e., essentially reaches its thermodynamic-limit () value. This implies two important conclusions. First, the slow, subdiffusive transport appears to be a genuine property of a long chain on the ergodic side of the MBL transition and not just a finite-size effect. Second, the estimated value of the critical disorder is , i.e., considerably larger than the values suggested by the exact-diagonalization analysis of small systems. This estimate is a lower bound determined by considering the lower end of our error estimate with two standard deviations, plus a small systematic error shown as triangle symbols in Fig. 2 ^{2}

In order to understand how representative the average value of the imbalance is, we have also studied its fluctuations, . The root-mean-square amplitude of the fluctuations of the imbalance, , has been found to be very slowly varying with time for long times. This allows us to consider the probability density function , where . In this window, is, to an excellent approximation, Gaussian-distributed (see Fig. 1b), with a standard deviation roughly proportional to . This implies self-averaging of the imbalance. A similar behavior is found for a non-interacting, Anderson-localized system. The wide distributions for smaller system sizes complicate the accurate determination of around the transition, so that considering larger system sizes is very beneficial also from this point of view.

### iii.2 Entanglement and Schmidt gap

The entropy of entanglement characterizes the spread of correlations in the system.
On the ergodic side, power laws have been predicted for the growth of entanglement close to the transition Luitz et al. (2016): .
Another quantity of interest is the so-called Schmidt gap , which is defined as the difference between the two largest values of the entanglement spectrum (see Ref. Laflorencie (2016) and references therein) for a bipartition of the system into subsystems and : .
Here the entanglement eigenvalues ^{3}

To explain the connection between the Schmidt gap and the more frequently used von Neumann entropy of entanglement, we recall that the latter can be expressed in terms of the entanglement eigenvalues as follows :

(4) |

Here the eigenvalues are ordered in a descending way by convention. The approximate equality holds as long as the entanglement remains relatively low. For our problem, we find that the effect of the cutoff by the bond dimension is negligible for from the perspective of (see Fig. 3).

Contrary to the entanglement entropy, the evaluation of the Schmidt gap requires knowledge of only the first two entanglement eigenvalues and hence is less sensitive to the value of the bond dimension. A further appealing property of this quantity is that, for a thermalized system, for , while for a localized system it is expected that the Schmidt gap remains finite in the long-time limit, , or at most decays logarithmically. In this sense, the behavior of the Schmidt gap allows one to distinguish between the localized and delocalized phases in the same way as for the imbalance. The entanglement entropy and Schmidt gap for and representative choices of are depicted in Figs. 3 and 4, respectively, which also show the distributions of these quantities close to the final time considered here. The shape of the distributions found for the entropy of entanglement is in qualitative agreement with a detailed study of such distributions for eigenstates of small systems using exact diagonalization Luitz (2016). That these distributions are converged with respect to system size reflects the slow spread of entanglement close to the transition (even on the ergodic side).

We have determined the Schmidt gap by using a bipartition in the middle of the chain. Interestingly, we find that the averaged Schmidt gap shows in the delocalized phase a power-law decay similar to that as . The corresponding power-law exponent is shown in Fig. 5. The results demonstrate a striking qualitative similarity to those for the imbalance exponent , Fig. 2. However, the found large- values of are somewhat above those for . The corresponding estimate for the critical disorder of the MBL transition obtained as a point where the Schmidt gap exponent vanishes within error bars is , i.e., even larger than found from results for the imbalance. For the entropy of entanglement , however, we clearly see the effect of the cutoff at imposed by the bond dimension (see Fig. 4a, left panel), where the power-law behavior is disturbed around . Considering that the transport properties have nonetheless converged with , this indicates that the TDVP can indeed provide reasonable results even if there are significant cutoff effects in terms of the entropy. Moreover, the Schmidt gap, which still shows a clean power-law behavior at , appears to be less sensitive to this numerical cutoff.

## Iv Machine learning

To further corroborate our analysis, we apply methods from machine learning Bishop (2006), which has emerged recently as a powerful tool to analyze localization phenomena Schindler et al. (2017); Mano and Ohtsuki (2017); Venderley et al. (2018); van Nieuwenburg et al. (2017a); Hsu et al. (2018); Zhang et al. (2018), to our data obtained using the TDVP. We use two algorithms: a partially supervised approach that has previously been employed in Ref. Schindler et al. (2017), and a fully unsupervised method based on the “learning by confusion” scheme introduced in Ref. van Nieuwenburg et al. (2017b). The combination of traditional numerical analysis and machine learning is mutually reinforcing: understanding localization though machine learning amounts to learning machine learning through localization.

### iv.1 Supervised classification algorithm

For the supervised learning approach, we train a feed-forward neural network to distinguish data at two extremes of our data set: (delocalized) and (presumed to be localized). We choose a single hidden layer network with a activation function for the hidden layer (of size ) and a activation function for the output layer (of size , corresponding to the two classes we are distinguishing). As input data, we use time-series , , evaluated at equidistant time-steps taken from the full TDVP time evolution at and . We choose a simple cross entropy error function with regularization, cf. Ref. Schindler et al. (2017).

After convergence of the training set error, we apply the trained network on a test set containing data at and , finding a classification accuracy of larger than . We subsequently apply the trained network to TDVP time-series at intermediate disorder strengths , to determine the average confidence with which data corresponding to a given disorder strength is classified as either belonging to the delocalized () or localized () class of time-series we trained with. Successful training implies that for and for , while at intermediate disorder strengths quantifies how similar the time evolution is to the extreme values and . The results are shown in Fig. 6.

This approach clearly indicates the formation of a plateau for , suggesting the presence of a many-body localized phase over a range of disorder strengths, but its thermal counterpart gives way to an extended crossover region. In contrast, no plateau or a less pronounced one is present in the data for system sizes and , respectively, showing the importance of considering larger systems.

These results can be interpreted in analogy with the analysis of the power-law exponents and : over a range of disorder strengths, the imbalance does not decay and the spin densities “look similar” to the case , down to the critical disorder. However, the algorithm can still pick up differences between a large value of the power law and a small one, leading to a non-zero even in the delocalized regime, so that the behavior of and is qualitatively similar (see Fig. 6, rightmost panel). The plateau, within error bars, starts at ^{4}

### iv.2 Unsupervised confusion algorithm

In Ref. van Nieuwenburg et al. (2017b) a scheme called “learning by confusion” for the unsupervised detection of phase transitions for data ordered along a one-dimensional parameter space (here: the strength of disorder ) has been proposed. For this, an arbitrary parameter is fixed and a feed-forward neural network is trained (in a supervised fashion) to distinguish data with from data with . The accuracy of classification that the trained network achieves when applied to a test set is evaluated. The process is repeated for different choices of . The resulting function has two global maxima at and , where the accuracy is trivially 1 because all the data can be classified as belonging to one phase. In Ref. van Nieuwenburg et al. (2017b) it was observed that an additional local maximum of occurs when equals the location of a phase transition , as it would be easiest for the network to classify the data into two sets for this choice of separation. Thus, in presence of a single phase transition as a function of , the curve is expected to take a W-shape.

We apply this algorithm, using the same neural network architecture and the same type of input data as in the supervised case, while dividing the training data by . Each time we determine the average accuracy as the percentage of correctly classified sets (with respect to the division ) in a test set (which does not coincide with the training data). We observe that the results depend strongly on the initial conditions for the training, i.e., on the random initial choice of weights and biases for the network. Even for input data from the largest system sizes , no consistent W-shape emerges in . This is due to substantial fluctuations of for fixed between training runs, in particular around the putative transition region.

The fluctuations of resulting from the initial conditions of the training are, however, not random. We claim that they carry information that can be used to locate the transition. We observe that over a large region of , they follow a bimodal distribution as shown in Fig. 7. We interpret the two branches as instances of trained networks that identify localized and ergodic features, respectively
^{5}

## V Summary and Outlook

The quench dynamics of large disordered spin chains in the Heisenberg model has been investigated by means of the time-dependent variational principle for matrix product states. We have studied the long-time behavior of the imbalance and found a regime, occurring in a broad range of parameters of the system, with slow, yet finite transport. In terms of averaged transport, we find that the average, typical and median imbalance all lead to the same power-law exponent . Our results imply that the ergodic regime extends (at least) up to disorder . We observe a substantial shift of the exponent (and, as a result, of an estimate for the MBL transition point ) when we go from relatively small systems (that can be exactly diagonalized) to large systems with and . On the other hand, we do not see any significant difference between the results for and , which indicates saturation of the exponent . This favors the conclusion that the subdiffusive behavior is a true long-time asymptotic behavior. These findings have been supported by the results for the dynamics of the Schmidt gap in the entanglement spectrum, which shows a very similar behavior.

We have complemented a conventional analysis of the data (with fitting the average values to power laws) by a machine-learning analysis of the whole time dependences of individual spins in various realizations of disorder. We have chosen two approaches, a supervised and a novel unsupervised algorithm. The latter is based on learning by confusion, but also exploits the stochastic nature of the learning process. This approach provides a more accurate way to determine the position of the transition, and could be useful for wider applications in determining phase transitions. The results from both the supervised and unsupervised machine-learning methods support our conclusions concerning the bound for the critical disorder strength . Thus, we can reinforce machine learning by combining it with a “traditional” analysis of numerical data, verifying that machine learning tools can be reliably applied.

Our results, showing a very slow transport on the ergodic side of the MBL transition () support the expectation that the system looks essentially localized at the transition point () Pal and Huse (2010). From this point of view, the MBL transition has much in common with the Anderson transition on random regular graphs (RRG) Tikhonov et al. (2016). The latter problem is viewed as a toy model of the MBL transition, even though this connection is less precise for short-range interaction models than for those with power-law interactions Tikhonov and Mirlin (2018). A slow dynamics in the RRG model was recently studied numerically in Ref. Biroli and Tarzia (2017). From the perspective of the machine-learning approaches applied in this work, the notion that the transition point is localized in character manifests itself in that it is located near the edge of the “plateau” characterizing the MBL phase, rather than at the midpoint between the two regimes. In terms of the confusion-based algorithm, the transition is located at the crossing between “delocalized” and “localized” branches.

A slow dynamics near the transition and a localized character of the critical point are also qualitatively consistent with the avalanche mechanism of the transition developed recently in Refs. De Roeck et al. (2016); Thiery et al. (2017a, b). This scenario predicts an extended weakly delocalized regime above the “nominal” threshold and implies that previous numerical results might be tainted by finite-size effects. Specifically, within this scenario, static “ergodic spots” (spatial regions with anomalously weak disorder) embedded in the nominally localized phase Gopalakrishnan et al. (2015) thermalize the rest of the system below the true MBL transition at . In small systems, such spots are typically not found and the chains appear to be localized for ; however, with increasing the probability of finding ergodic spots increases and the true transition at becomes well resolved.

Indeed, we find that considering larger systems leads to a substantial increase in delocalization, in agreement with the avalanche scenario. Moreover, the unsupervised machine-learning analysis, showing the coexistence of localizing and delocalizing realizations of disorder in the range (see Fig. 7), may also be interpreted in terms of the avalanche-induced delocalization. Further, consistent with the prediction of Thiery et al. Thiery et al. (2017a, b), the system “looks localized” near the transition, as evidenced by the strong similarity between the spin dynamics for the non-interacting case and the interacting case close to our lower bound for the MBL transition. However, at the present stage, we cannot explicitly confirm (or falsify) this mechanism, as this would require a specific analysis of observables that would have a distinct behavior within this mechanism. Another possibility consistent with our results is the existence of a glass-like crossover regime Biroli and Tarzia (2017).

Future work inspired by our results could focus on different energy densities in addition to just the center of the band that has been considered here, in order to map the transition line in the energy-disorder plane on the basis of data for large systems. Further, it is also interesting and instructive to analyze the phase diagram of long XXZ chains in the interaction-disorder plane. Clearly, a better analytical understanding of the slow dynamics near the MBL transition would be very important. In this context, it is interesting to note that a similar slow transport near the transition (or, more accurately, apparent transition) is also found for quasiperiodic systems Schreiber et al. (2015) as well as in two-dimensional disordered systems Bordia et al. (2017), where the influence of rare bottlenecks is expected to be negligible. It remains to be seen whether this slow transport has a common origin in all these situations.

A so-called “dreaming” protocol within the neural-network framework, which generates the configurations that are representative for given phases, could be very useful for developing the analytical theories of the subdiffusive phase. Furthermore, this type of machine-learning approach can be envisaged to simulate the dynamics at much longer times (and, perhaps, in much larger systems), which is currently inaccessible by the most advanced numerical tools. Finally, a hotly debated question is the very existence of the true () MBL transition in various models. We hope that our work will pave the way for numerical investigations of the physics associated with the MBL transition on large systems, which is of obvious importance for shedding light on this issue.

###### Acknowledgements.

We thank F. Alet, Y. Bar Lev, M. H. Fischer, C. Karrasch, N. Laflorencie, D. Luitz, S. R. Manmana, M. Müller, and R. M. Nandkishore for useful discussions. TDVP simulations were performed using a Python script based on the open-source evoMPS library Milsted et al. (2013), implementing the single-site algorithm of Ref. Haegeman et al. (2016). We are grateful to A. Milsted for help with implementing evoMPS. TEBD simulations were performed using Open Source Matrix Product States Wall and Carr (2012); Jaschke et al. (2018). Simulations for the non-interacting model were performed using the NumPy implementation of LAPACK Anderson et al. (1999). Machine learning algorithms were implemented using TensorFlow Abadi et al. (2015). We used Matplotlib Hunter (2007) to generate figures and GNU Parallel Tange (2011) for running parallel simulations. The authors acknowledge support by the state of Baden-Württemberg through bwHPC. The work of KST, IVG, and ADM was supported by the Russian Science Foundation (Grant No. 14-42-00044). FS and TN acknowledge support from the Swiss National Science Foundation (grant number: 200021_169061) and from the European Union’s Horizon 2020 research and innovation program (ERC-StG-Neupert-757867-PARATOP).## Appendix A Convergence of the TDVP

Since the TDVP truncates the entanglement spectrum, it is of interest to investigate the convergence with respect to bond dimension Kloss et al. (2018). We perform two benchmarks: first, we compare the time evolution of the TDVP to numerically exact numerics for a small system . Here, we consider the disorder average of independent disorder samples, since we are interested in convergence at the level of this average, rather than at the level of individual realizations. The exact numerics are obtained simply by performing TDVP time evolution with an unrestricted bond dimension. For a system of size an MPS with maximum bond dimension captures the time evolution exactly (for a system of size , non-truncated MPS have a maximum bond dimension of in the center of the chain). The result is shown in Fig. 8, which shows the time evolution of the initial state up to times . Recall that in the main body of the paper, we considered times only up to . Secondly, we consider a larger system and consider whether the result converges with bond dimension, where we consider times up to and bond dimensions up to , higher than in the main body of the paper. The result is shown in Fig. 9.

From both of these approaches, we conclude that for the times we consider (up to ), a bond dimension of provides an excellent approximation for not too weak disorder ; an approximation that only improves with increasing disorder. This conclusion is further corroborated by comparison with a completely independent implementation of the time-evolving block decimation method, for which we compute the power law decay coefficient for , finding excellent agreement with the TDVP (see Fig. 2). Therefore, our numerical results are certainly accurate in the disorder and time window we have considered in this work. Moreover, in the moderately strongly disordered case we find good agreement between power laws obtained in the window and those in the window , providing evidence for the survival of these power laws in the long-time limit.

### Footnotes

- It should be noted that for small , we can no longer distinguish a power-law decay over our limited time window from, say, a logarithmic correction. Nevertheless, this procedure yields a sensitive probe for the decay of in the sense that any statistically significant monotonous decay will result in a non-zero .
- The imbalance of a non-interacting system , while asymptotically approaching a positive constant, shows a weak decay with time. We have checked that fitting of this decay by a power law in the considered time interval yields fictitious values of in the range . This implies a lower bound on the true exponent of a many-body delocalized system that we can distinguish from zero.
- In numerical language, the entanglement eigenvalues are also the squares of the Schmidt numbers, obtained easily in the MPS representation through singular value decompositions.
- The start of the plateau can be associated with the value of where the lower bound of the error bar in the right panel of Fig. 6 touches the upper bound of the error bar for the training at .
- It is interesting to note that the “MBL branch” for large chains ( and ) starts at , which is close to the estimate for the critical disorder that was previously obtained from the exact diagonalization in relatively small chains.

### References

- Rahul Nandkishore and David A. Huse, “Many-body localization and thermalization in quantum statistical mechanics,” Ann. Rev. Cond. Mat. Phys. 6, 15–38 (2015).
- Ehud Altman and Ronen Vosk, “Universal dynamics and renormalization in many-body-localized systems,” Ann. Rev. Cond. Mat. Phys. 6, 383–409 (2015).
- Dmitry A. Abanin and Zlatko Papić, “Recent progress in many-body localization,” Ann. Phys. (Berl.) 529, 1700169 (2017).
- Ferdinand Evers and Alexander D. Mirlin, “Anderson transitions,” Rev. Mod. Phys. 80, 1355 (2008).
- M. Rigol, V. Dunjko, and M. Olshanii, “Thermalization and its mechanism for generic isolated quantum systems,” Nature 452, 854 (2008).
- I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, “Interacting electrons in disordered wires: Anderson localization and low- transport,” Phys. Rev. Lett. 95, 206603 (2005).
- D.M. Basko, I.L. Aleiner, and B.L. Altshuler, “Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states,” Ann. Phys. (N. Y.) 321, 1126 – 1205 (2006).
- Marko Žnidarič, Antonello Scardicchio, and Vipin Kerala Varma, “Diffusive and subdiffusive spin transport in the ergodic phase of a many-body localizable system,” Phys. Rev. Lett. 117, 040601 (2016).
- Vedika Khemani, Frank Pollmann, and S. L. Sondhi, “Obtaining highly excited eigenstates of many-body localized hamiltonians by the density matrix renormalization group approach,” Phys. Rev. Lett. 116, 247204 (2016).
- Xiongjie Yu, David Pekker, and Bryan K. Clark, “Finding matrix product state representations of highly excited eigenstates of many-body localized hamiltonians,” Phys. Rev. Lett. 118, 017201 (2017).
- Thorsten B. Wahl, Arijeet Pal, and Steven H. Simon, “Efficient representation of fully many-body localized systems using tensor networks,” Phys. Rev. X 7, 021018 (2017).
- Elmer V. H. Doggen, Gabriel Lemarié, Sylvain Capponi, and Nicolas Laflorencie, “Weak- versus strong-disorder superfluid—bose glass transition in one dimension,” Phys. Rev. B 96, 180202 (2017).
- Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, Henrik P. Lüschen, Mark H. Fischer, Ronen Vosk, Ehud Altman, Ulrich Schneider, and Immanuel Bloch, “Observation of many-body localization of interacting fermions in a quasirandom optical lattice,” Science 349, 842–845 (2015).
- Jae-yoon Choi, Sebastian Hild, Johannes Zeiher, Peter Schauß, Antonio Rubio-Abadal, Tarik Yefsah, Vedika Khemani, David A. Huse, Immanuel Bloch, and Christian Gross, “Exploring the many-body localization transition in two dimensions,” Science 352, 1547–1552 (2016).
- J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, ‘‘Many-body localization in a quantum simulator with programmable random disorder,” Nat. Phys. 12, 907 (2016).
- G. Kucsko, S. Choi, J. Choi, P. C. Maurer, H. Zhou, R. Landig, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, E. Demler, N. Y. Yao, and M. D. Lukin, arXiv:1609.08216 .
- John Z. Imbrie, “On many-body localization for quantum spin chains,” J. Stat. Phys. 163, 998–1048 (2016).
- Arijeet Pal and David A. Huse, “Many-body localization phase transition,” Phys. Rev. B 82, 174411 (2010).
- David J. Luitz, Nicolas Laflorencie, and Fabien Alet, ‘‘Many-body localization edge in the random-field heisenberg chain,” Phys. Rev. B 91, 081103 (2015).
- F. Pietracaprina, N. Macé, D. J. Luitz, and F. Alet, arXiv:1803.05395 (2018).
- Timothy C. Berkelbach and David R. Reichman, “Conductivity of disordered quantum lattice models at infinite temperature: Many-body localization,” Phys. Rev. B 81, 224429 (2010).
- Trithep Devakul and Rajiv R. P. Singh, “Early breakdown of area-law entanglement at the many-body delocalization transition,” Phys. Rev. Lett. 115, 187201 (2015).
- Maksym Serbyn, Z. Papić, and Dmitry A. Abanin, “Universal slow growth of entanglement in interacting strongly disordered systems,” Phys. Rev. Lett. 110, 260601 (2013).
- Yevgeny Bar Lev, Guy Cohen, and David R. Reichman, “Absence of diffusion in an interacting system of spinless fermions on a one-dimensional disordered lattice,” Phys. Rev. Lett. 114, 100601 (2015).
- Kartiek Agarwal, Sarang Gopalakrishnan, Michael Knap, Markus Müller, and Eugene Demler, “Anomalous diffusion and griffiths effects near the many-body localization transition,” Phys. Rev. Lett. 114, 160401 (2015).
- Kartiek Agarwal, Ehud Altman, Eugene Demler, Sarang Gopalakrishnan, David A. Huse, and Michael Knap, “Rare-region effects and dynamics near the many-body localization transition,” Ann. Phys. (Berl.) 529, 1600326 (2016).
- David J. Luitz, Nicolas Laflorencie, and Fabien Alet, “Extended slow dynamical regime close to the many-body localization transition,” Phys. Rev. B 93, 060201 (2016).
- Soumya Bera, Giuseppe De Tomasi, Felix Weiner, and Ferdinand Evers, “Density propagator for many-body localization: Finite-size effects, transient subdiffusion, and exponential decay,” Phys. Rev. Lett. 118, 196801 (2017).
- P. A. M. Dirac, “Note on exchange phenomena in the thomas atom,” Proc. Cambridge Philos. Soc. 26, 376–385 (1930).
- Jutho Haegeman, J. Ignacio Cirac, Tobias J. Osborne, Iztok Pižorn, Henri Verschelde, and Frank Verstraete, “Time-dependent variational principle for quantum lattices,” Phys. Rev. Lett. 107, 070601 (2011).
- Jutho Haegeman, Christian Lubich, Ivan Oseledets, Bart Vandereycken, and Frank Verstraete, “Unifying time evolution and optimization with matrix product states,” Phys. Rev. B 94, 165116 (2016).
- A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, “Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces,” J. Stat. Mech.: Th. Exp. 2004, P04005 (2004).
- Steven R. White and Adrian E. Feiguin, ‘‘Real-time evolution using the density matrix renormalization group,” Phys. Rev. Lett. 93, 076401 (2004).
- G. Vidal, “Efficient classical simulation of slightly entangled quantum computations,” Phys. Rev. Lett. 91, 147902 (2003).
- Pasquale Calabrese and John Cardy, “Evolution of entanglement entropy in one-dimensional systems,” J. Stat. Mech.: Th. Exp. 2005, P04010 (2005).
- Nicolas Laflorencie, Phys. Rep. 646, 1 – 59 (2016).
- U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Ann. Phys. (N. Y.) 326, 96 – 192 (2011).
- C. Karrasch, J. H. Bardarson, and J. E. Moore, “Reducing the numerical effort of finite-temperature density matrix renormalization group calculations,” New J. Phys. 15, 083031 (2013).
- E. Leviatan, F. Pollmann, J. H. Bardarson, D. A. Huse, and E. Altman, arXiv:1702.08894 (2017).
- Benedikt Kloss, Yevgeny Bar Lev, and David Reichman, “Time-dependent variational principle in matrix-product state manifolds: Pitfalls and potential,” Phys. Rev. B 97, 024307 (2018).
- Mario Berta, Fernando G. S. L. Brandão, Jutho Haegeman, Volkher B. Scholz, and Frank Verstraete, arXiv:1709.07423 (2017).
- J. M. Deutsch, “Quantum statistical mechanics in a closed system,” Phys. Rev. A 43, 2046–2049 (1991).
- M. Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
- David J. Luitz and Yevgeny Bar Lev, “The ergodic side of the many-body localization transition,” Ann. Phys. (Berl.) 529, 1600350–n/a (2017).
- Daniel Jaschke, Michael L. Wall, and Lincoln D. Carr, “Open source matrix product states: Opening ways to simulate entangled many-body quantum systems in one dimension,” Comput. Phys. Commun. 225, 59 – 91 (2018).
- It should be noted that for small , we can no longer distinguish a power-law decay over our limited time window from, say, a logarithmic correction. Nevertheless, this procedure yields a sensitive probe for the decay of in the sense that any statistically significant monotonous decay will result in a non-zero .
- The imbalance of a non-interacting system , while asymptotically approaching a positive constant, shows a weak decay with time. We have checked that fitting of this decay by a power law in the considered time interval yields fictitious values of in the range . This implies a lower bound on the true exponent of a many-body delocalized system that we can distinguish from zero.
- In numerical language, the entanglement eigenvalues are also the squares of the Schmidt numbers, obtained easily in the MPS representation through singular value decompositions.
- Hui Li and F. D. M. Haldane, “Entanglement spectrum as a generalization of entanglement entropy: Identification of topological order in non-abelian fractional quantum hall effect states,” Phys. Rev. Lett. 101, 010504 (2008).
- David J. Luitz, “Long tail distributions near the many-body localization transition,” Phys. Rev. B 93, 134201 (2016).
- Christopher Bishop, Pattern Recognition and Machine Learning (Springer-Verlag, New York, 2006).
- Frank Schindler, Nicolas Regnault, and Titus Neupert, “Probing many-body localization with neural networks,” Phys. Rev. B 95, 245134 (2017).
- Tomohiro Mano and Tomi Ohtsuki, “Phase diagrams of three-dimensional anderson and quantum percolation models using deep three-dimensional convolutional neural network,” J. Phys. Soc. Jpn. 86, 113704 (2017).
- Jordan Venderley, Vedika Khemani, and Eun-Ah Kim, “Machine learning out-of-equilibrium phases of matter,” Phys. Rev. Lett. 120, 257204 (2018).
- E. P. L. van Nieuwenburg, E. Bairey, and G. Refael, arXiv:1712.00450 (2017a).
- Y.-T. Hsu, X. Li, D.-L. Deng, and S. Das Sarma, arXiv:1805.12138 (2018).
- Wei Zhang, Lei Wang, and Ziqiang Wang, arXiv:1807.02954 (2018).
- Evert P. L. van Nieuwenburg, Ye-Hua Liu, and Sebastian D. Huber, “Learning phase transitions by confusion,” Nat. Phys. 13, 435 (2017b).
- The start of the plateau can be associated with the value of where the lower bound of the error bar in the right panel of Fig. 6 touches the upper bound of the error bar for the training at .
- It is interesting to note that the “MBL branch” for large chains ( and ) starts at , which is close to the estimate for the critical disorder that was previously obtained from the exact diagonalization in relatively small chains.
- K. S. Tikhonov, A. D. Mirlin, and M. A. Skvortsov, “Anderson localization and ergodicity on random regular graphs,” Phys. Rev. B 94, 220203 (2016).
- K. S. Tikhonov and A. D. Mirlin, “Many-body localization transition with power-law interactions: Statistics of eigenstates,” Phys. Rev. B 97, 214205 (2018).
- G. Biroli and M. Tarzia, “Delocalized glassy dynamics and many-body localization,” Phys. Rev. B 96, 201114 (2017).
- Wojciech De Roeck, Francois Huveneers, Markus Müller, and Mauro Schiulaz, ‘‘Absence of many-body mobility edges,” Phys. Rev. B 93, 014203 (2016).
- T. Thiery, F. Huveneers, M. Müller, and W. De Roeck, arXiv:1706.09338 (2017a).
- T. Thiery, M. Müller, and W. De Roeck, arXiv:1711.09880 (2017b).
- Sarang Gopalakrishnan, Markus Müller, Vedika Khemani, Michael Knap, Eugene Demler, and David A. Huse, “Low-frequency conductivity in many-body localized systems,” Phys. Rev. B 92, 104202 (2015).
- Pranjal Bordia, Henrik Lüschen, Sebastian Scherg, Sarang Gopalakrishnan, Michael Knap, Ulrich Schneider, and Immanuel Bloch, “Probing slow relaxation and many-body localization in two-dimensional quasiperiodic systems,” Phys. Rev. X 7, 041047 (2017).
- Ashley Milsted, Jutho Haegeman, Tobias J. Osborne, and Frank Verstraete, “Variational matrix product ansatz for nonuniform dynamics in the thermodynamic limit,” Phys. Rev. B 88, 155116 (2013).
- M. L. Wall and Lincoln D. Carr, “Out-of-equilibrium dynamics with matrix product states,” New J. Phys. 14, 125015 (2012).
- E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users’ Guide (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
- Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,” (2015), software available from tensorflow.org.
- J. D. Hunter, “Matplotlib: A 2d graphics environment,” Comput. Sci. Eng. 9, 90–95 (2007).
- O. Tange, “Gnu parallel - the command-line power tool,” ;login: The USENIX Magazine 36, 42–47 (2011).