{'fraction_non_alphanumeric': 0.044385609718112444, 'fraction_numerical': 0.017442765924310856, 'mean_word_length': 4.529276693455798, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Domain-wall melting as a probe of many-body localization

Motivated by a recent optical-lattice experiment by Choi et al.[Science 352, 1547 (2016)], we discuss how domain-wall melting can be used to investigate many-body localization. First, by considering noninteracting fermion models, we demonstrate that experimentally accessible measures are sensitive to localization and can thus be used to detect the delocalization-localization transition, including divergences of characteristic length scales. Second, using extensive time-dependent density matrix renormalization group simulations, we study fermions with repulsive interactions on a chain and a two-leg ladder. The extracted critical disorder strengths agree well with the ones found in existing literature.

Introduction. In pioneering works based on perturbation theory [1,2], it was shown that Anderson localization, i.e., perfectly insulating behavior even at finite temperatures, can persist in the presence of interactions. Subsequent theoretical studies on mostly onedimensional (1D) model systems have unveiled many fascinating properties of such a many-body localized (MBL) phase. The MBL phase is a dynamical phase of matter defined in terms of the properties of highly excited manybody eigenstates. It is characterized by an area-law entanglement scaling in all eigenstates [3][4][5], a logarithmic increase of entanglement in global quantum quenches [6][7][8], failure of the eigenstate thermalization hypothesis [9] and therefore, memory of initial conditions [10,11]. The phenomenology of MBL systems is connected to the existence of a complete set of commuting (quasi) local integrals of motion (so-called "l-bits") that are believed to exist in systems in which all many-body eigenstates are localized [8,12,13]. These l-bits can be thought of as quasiparticles with an infinite lifetime, in close analogy to a zero-temperature Fermi liquid [1,14]. Important open questions pertain to the nature of the MBL transition and the existence of an MBL phase in higher dimensions, for which there are only few results (see, e.g., [15,16]), mainly due to the fact that numerical simulations are extremely challenging in dimensions higher than one for the MBL problem.
The phenomenology of the MBL phase has mostly been established for closed quantum systems. A sufficiently strong coupling of a disordered, interacting system to a bath is expected to lead to thermalization (see, e.g., [17,18]). Thus, the most promising candidate systems for the experimental investigation of MBL physics are quantum simulators such as ultracold quantum gases in optical lattices or ion traps. So far, the cleanest evidence for MBL in an experiment has been reported for an interacting Fermi gas in an optical lattice with quasiperiodicity, realizing the Aubry-André model [19,20]. Other quantum gas experiments used the same quasi-periodic lattices or laser speckles to investigate Anderson localization [21,22] and the effect of interactions [23], however, at low energy densities. Experiments with ion traps provide an alternative route, yet there, at most a dozen of ions can currently be studied [24].
By using a novel experimental approach, a first demonstration and characterization of MBL in a twodimensional (2D) optical-lattice system of interacting bosons with disorder has been presented by Choi et al. [25]. They start from a state that contains particles in only one half of the system while the rest is empty. Once tunneling is allowed, the particles from the initially occupied region can spread out into the empty region (see Fig. 1). The evolution of the particle density is tracked using single-site resolution techniques [26,27] and digital mirror devices allow one to tune the disorder. The relaxation dynamics provides evidence for the existence of an ergodic and an MBL regime as disorder strength is varied, characterized via several observables such as density profiles, particle-number imbalances and measures of the localization length [25]. This experiment serves as the main motivation for our theoretical work.
The term domain-wall melting is inherited from the equivalent problem in quantum magnetism (see, e.g., [28][29][30][31][32][33]), corresponding to coupling two ferromagnetic domains with opposite spin orientation. Furthermore, the domain-wall melting describes the transient dynamics [34][35][36] of sudden-expansion experiments of interacting quantum gases in optical lattices (i.e., the release of initially trapped particles into an empty homogeneous lattice) [36][37][38][39]. Theoretically, the sudden expansion of interacting bosons in the presence of disorder was stud- ied in, e.g., [40,41] for the expansion from the correlated ground state in the trap, while for MBL, higher energy densities are relevant.
We use exact diagonalization (ED) and timedependent density matrix renormalization group (tDMRG) methods [42][43][44][45] to clarify some key questions of the domain-wall experiments. First, by considering noninteracting fermions in a 1D tight-binding model with diagonal disorder we demonstrate that it is possible to extract the single-particle localization length ξ (1) loc as a function of disorder strength from such an experiment since the density profiles develop exponential tails with a length scale ξ dw (see Fig. 1). This domain-wall decay length ξ dw also captures the disorder driven metal-insulator transition in the Aubry-André model when approached from the localized regime, exhibiting a divergence. Second, we study the case of spinless fermions with nearest-neighbor repulsive interactions on chains and two-leg ladders, for which numerical estimates of the critical disorder strength W c of the metal-insulator transition are available [5,9,14,[46][47][48][49]. For both models, essential features of the noninteracting case carry over, namely, the steady-state profiles decay exponentially with distance in the localized regime W > W c (i.e., the expansion stops), while particles continue to spread in the ergodic regime W < W c . Moreover, we discuss experimentally accessible measures to investigate the dynamics close to the transition for all models.
Noninteracting cases. We start by considering fermions in a 1D tight-binding lattice with uncorrelated diagonal disorder. The Hamiltonian reads: whereĉ † j denotes the creation operator on site j,n j = c † jĉj is the number operator, n j = n j is density, and j ∈ [−W, W ] is a random onsite potential (L is the number of sites). We set the lattice spacing to unity and = 1. All single-particle eigenstates are localized for any nonzero W and thus the system is an Anderson insulator at all energy densities [50,51]. Typical density profiles for the dynamics starting from a domain-wall initial state are shown for different times in Fig. 2(a). Here, "typical" refers to the geometric mean n j over disorder realizations (i.e., the arithmetic mean of log n j ) (see [52]). The domain wall first melts slightly yet ultimately stops expanding. The profiles clearly develop an exponential tailn j ∝ exp(−2j/ξ dw ) for j 0. The crucial question is now whether the length scale ξ dw is directly related to the single-particle localization length or not.
We compare two ways of extracting ξ dw : First, a fit to the numerical data forn j in the tails j 0 and second, via computing the variance of the particles emit-  ted into the originally empty region. For the latter, we view the density n j in the initially empty region j > 0 as a spatial distribution · n ≡ j>0 n j · /∆N where ∆N = j>0 n j is the number of emitted particles. The variance VAR n = j 2 n − j 2 n of this particle distribution is shown in Fig. 2(c) and approaches a stationary regime on a timescale depending on W . For the time window plotted, only the curves with W ≥ J saturate, yet we checked that also the curves for W < J saturate at sufficiently long times. At short times, VAR n ∝ t 2 signals a ballistic expansion of the particles as long as VAR n (t) ξ loc . Assuming a strictly exponential distribution n j ∝ exp(− 2j ξ dw ) for all j > 0 yields VAR n ≈ ξ 2 dw 4 for VAR n 1. We use that relation to extract ξ dw in the general case as well and in addition, we introduce an explicit time dependence of ξ dw to illustrate the approach to the stationary state. In general, this gives only a lower bound to W c since VAR n can be finite for diverging ξ dw if the distribution is not exponential. Yet we find that both methods give similar results for the final profile and show only ξ dw extracted from VAR n in Fig. 2(b).
The known result for the localization length in the 1D Anderson model is ξ [50] for E = 0 (our initial state leads to that average energy for sufficiently large systems). Our data for ξ dw shown in Fig. 2(b) clearly exhibit the expected scaling ξ dw ∝ W −2 over a wide range of W as suggested by a fit of ξ dw = a/W −2 to the data [dashed line in Fig. 2(b); the prefactor is larger by about a factor of 1.5 than the typical localization length ξ (1) loc ]. Deviations from the W −2 dependence at small W , where ξ (1) loc ∼ O(L), are due to the finite system size. At large W , the discreteness of the lattice makes it impossible to resolve ξ dw that are much smaller than the lattice spacing. We stress that fairly long times need to be reached to observe a good quantitative agreement with the W −2 dependence. For instance, for the parameters of Fig. 2(a), tJ ∼ 1000 is necessary to reach the asymptotic form. Nevertheless, even at shorter times, the density profiles are already approximately exponential. To summarize, our results demonstrate that the characteristic length scale ξ dw is a measure of the singleparticle localization length, most importantly exhibiting the same qualitative behavior.
In Fig. 2(d), we introduce an alternative indicator of localization, namely, the number of emitted particles ∆N (t) that have propagated across the edge j = 0 of the initial domain wall at a time t. Due to particle conservation, ∆N is directly related to the imbalance I = N −2∆N N analyzed in the experiment [25]. We observe that ∆N shares qualitatively the same behavior with VAR n [note the linear y scale in Fig. 2(d)], which will also apply to the models discussed in the following.
As a further test, we study the Aubry-André model in the Appendix [52]. The comparison of ξ dw with the exactly known single-particle localization length [53] in the inset of Fig. 2(b) demonstrates that the domain-wall melting can resolve the delocalization-localization transition at W = J.
Interacting fermions on a chain. Given the encouraging results discussed above, we move on to studying the dynamics in a system with an MBL phase, namely to the model of spinless fermions with repulsive nearestneighbor interactions H int = H 0 + V i,j n inj , equivalent to the spin-1/2 XXZ chain. We focus on SU (2) symmetric exchange, i.e., V = J, for which numerical studies predict a delocalization-localization transition from an ergodic to the MBL phase at W c /J = (3.5 ± 1) [5,9,14,47] at energy densities in the middle of the many-body spectrum (corresponding to infinite temperature when approaching the transition from the ergodic side). Note, though, that even for this much studied model, some aspects of the phase diagram are still de- bated in the recent literature (see, e.g., [54,55]). Typical time evolutions of density profiles in the ergodic and MBL phase are shown in Figs. 3(a)-3(c), obtained from tDMRG simulations [42][43][44]. We use a time step of dt = 0.04/J and a bond dimension of up to χ = 1000 and keep the discarded weight in each time step under 10 −10 . The disorder average is performed over about 500 realizations. These profiles show a crucial difference between the dynamics in the localized and the delocalized regime. Deep in the localized regime, Fig. 3(c), similar to the noninteracting models discussed before, the density profiles quickly become stationary with an exponential decay even close to j = 0. In the ergodic phase, however, the density profiles never become stationary on the simulated time scales and for the values of interactions considered here. For W = 0.5J shown in Fig. 3(a), the particles spread over the whole considered system. Remarkably, we find a regime of slow dynamics [47,[56][57][58][59] at intermediate disorder W < W c in Fig. 3(b), where there seems to persist an exponential decay of n j at finite times, but with a continuously growing ξ dw (t). We note that ξ dw (t) at the shortest time scales is on the order of the single-particle localization length. An explana- tion can thus be obtained in this picture: On short time scales, single particles can quickly expand into the right, empty side within the single-particle localization length, thus leading to the exponential form of n j . The interaction comes into play by scattering events at larger times, ultimately allowing the expansion over the whole system for infinite times.
The slow regime is also reflected in the quantities VAR n and ∆N in Figs. 3(d) and 3(e), which behave qualitatively in the same way. While both quantities saturate for W > W c and the results hardly differ from the noninteracting case shown by the dotted lines, the slow growth becomes evident for W W c at the intermediate time scales accessible to us. The slow growth of both VAR n and ∆N is, for W W c , the best described by (yet hard to distinguish from a power-law) ∆N (t), VAR n (t) = C(W ) log(tJ) + const . ( This growth is qualitatively different from the noninteracting case, where a saturation sets in after a faster initial increase. The inset of Fig. 3(e) shows the prefactor C(W ) extracted from a fit to the data of ∆N (t) for tJ > 10. This allows us to extract W c since C(W > W c ) = 0 for the stationary profiles in the localized phase. Our result for W c is compatible with the literature value W c /J = 3.5 ± 1 [5,9,14,47] (dashed line in the inset of Fig. 3).
Interacting fermions on a ladder. As a first step towards 2D systems, we present results for the dynamics of interacting spinless fermions on a two-leg ladder in the presence of diagonal disorder. The simulations are done with a variant of tDMRG suitable for long-range interactions [45], with a time step dt = 0.01/J. Figures  4(a) and 4(b) show the variance VAR n and ∆N for V = J, respectively. As for the chain, we observe that both the variance and ∆N have a tendency to saturate for large disorder strength, while they keep growing for small disorder. The data are best described by Eq. (2) and we extract C(W ) from fits of the data for tJ > 10 to Eq. (2). The results of these fits shown in the inset of Fig. 4 suggest a critical disorder strength 8 W c /J 10, in good agreement with the value of W c /J = 8.5 ± 0.5 found in an ED study of the isotropic Heisenberg model on a two-leg ladder [49] (the two models differ by correlated hopping terms which are not believed to be important for the locus of the transition).
Summary and outlook. We analyzed the domainwall melting of fermions in the presence of diagonal disorder, motivated by a recent experiment [25] that was first in using this setup for interacting bosons in 2D. Our main result is that several quantities accessible to experimentalists (such as the number of propagating particles and the variance of their particle density) are sensitive to localization and can be used to locate the disorder-driven metal-insulator transition, based on our analysis of several models of noninteracting and interacting fermions for which the phase diagrams are known. Notably, this encompasses a two-leg ladder as a first step towards numerically simulating the dynamics of interacting systems with disorder in the 1D-2D crossover. Our work further indicates that care must be taken in extracting quantitative results from finite systems or finite times since the approach to the stationary regime can be slow. Interestingly, we observe a slow dynamics in the ergodic phase of interacting models as the transition to the MBL phase is approached, which deserves further investigation.
The domain-wall melting thus is a viable approach for theoretically and experimentally studying disordered interacting systems, and we hope that our work will influence future experiments on quasi-1D systems where a direct comparison with theory is feasible. Concerning 2D systems, where numerical simulations of real-time dynamics face severe limitations, our results for two-leg ladders provide confidence that the domain-wall melting is still a reliable detector of localization as well, as evidenced in the experiment of [25]. Even for clean systems, experimental studies of domain-wall melting in the presence of interactions could provide valuable insights into the nonequilibrium transport properties of interacting quantum gases [28,29,34,35,37,38]. For instance, even for the isotropic spin-1/2 chain (V = 1 in our case), the qualitative nature of transport is still an open issue [60][61][62][63][64][65][66][67]. Moreover, the measurement of diffusion constants would be desirable [68].

Disorder Statistics
An exemplary distribution of n j for the free fermion case is shown in Fig. S1. In a rough approximation, the probability for a particle to hop the j sites out of the domain wall can be seen as a product of the hopping probabilities to neighboring sites, which depend on the specific disorder realization. The geometric meann j is thus a natural choice for the average over different disorder realizations. As evident from Fig. S1, it coincides with the median and represents the typical value. In contrast, the arithmetic mean is an order of magnitude larger as it puts a large weight in the upper tail of the distribution.
Although the geometric mean is a good choice for n j , it is reasonable to use the arithmetic mean for other quantities such as VAR n and ∆N : they represent quantities integrated over j for a given disorder realization. We checked that the arithmetic mean is close to typical values for these quantities.

Aubry-André model
We now focus on the dynamics in the Aubry-André model, where a quasiperiodic modulation is introduced in Eq. (1) via j = W cos(2πrj + φ 0 ) (employed in the MBL experiments of [19,20]). We set the irrational ratio r to r = ( √ 5−1)/2 = 0.61803 . . . and perform the equivalent to disorder averages by sampling over the value of the phase φ 0 ∈ [0, 2π). This noninteracting model has a delocalization-localization transition at W c /J = 1, where the single-particle localization length diverges as ξ (1) loc = log W J [53]. Similar to the previously considered Anderson model, the density profiles become stationary with an exponential tail in the localized phase for W > W c . As W is varied, a clear transition is visible in the time dependence of both VAR n and ∆N shown in Figs. S2(a) and S2(b), respectively, which become stationary for W > W c , while growing with a power law for W < W c . The corresponding domain-wall decay length ξ dw (see the inset of Fig. 2(b)) diverges as W c is approached from above, in excellent agreement with the single-particle localization length of that model [53]. The maximum value of ξ dw in the extended phase reached at long times diverges with L.